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19.  ABSTRACT  (CONCLUDED) 

cases  with  a  shock  wave  passing  through  the  block  boundary.  Converged  steady  results 
with  no  body  motion,  show  no  error  due  to  the  block  boundaries.  Convergence  is 
slowed  to  some  extent  by  using  blocked  grids,  but  the  use  of  blocked  grids  is 
economically  attractive  due  to  the  greatly  reduced  computer  memory  cost.  Steady  and 
unsteady  calculations  of  transonic  flow  over  a  rectangular  planform  supercritical 
wing  using  a  blocked  grid  are  compared  to  experimental  results.  The  comparison  is 
good  for  those  cases  in  which  viscous  effects  are  negligible,  but  for  higher  Mach 
number  cases  the  inviscid  calculation  places  the  shock  downstream  of  the  experimental 
location. 
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CHAPTER  I 


s 


INTRODUCTION 


Accurate  prediction  of  unsteady  transonic  flow  about  modern  fighter 


aircraft  is  important  due  to  the  requirement  for  transonic  combat  and 
store  separation,  and  because  of  the  aerodynamic  and  aeroelastic  prob¬ 
lems  encountered  in  this  regime.  Flutter,  for  example,  is  a  potentially 
destructive  coupling  between  flexible  structural  motion  and  the  unsteady 
aerodynamics  about  that  structure  that  is  typically  more  critical  in  the 
transonic  flight  regime  than  in  either  the  subsonic  or  supersonic.  Also, 
the  trajectory  of  a  store  separating  from  the  parent  aircraft  can  be 
greatly  affected  by  the  increased  aerodynamic  interference  between  air¬ 
craft  and  store  that  i3  typical  of  the  transonic  regime.  Although  these 
problems  have  been  recognized  for  some  time,  progress  in  the  prediction 
of  unsteady  transonic  aerodynamics  about  complex  geometries,  such  as  a 
fighter  aircraft  with  external  stores,  has  been  slow  due  to  the  inherent 
nonlinearity  of  transonic  flow  and  the  difficulties  involved  in  obtain¬ 
ing  solutions  to  the  equations  describing  this  flow. 

The  Navier-Stokes  equations  adequately  describe  unsteady  transonic 
flow,  including  shock  waves,  viscous  effects,  and  shock  and  boundary 
layer  interactions.  Solution  of  the  Navier-Stokes  equations  about  com¬ 
plex  three-dimensional  geometries  on  present-day  computers  would  be  very 
expensive.  Fortunately,  many  flow-fields  of  practical  interest  are  mod¬ 
eled  well  by  the  Euler  equations,  which  may  be  obtained  from  the 
Navier-Stokes  equations  by  assuming  that  the  flow  is  inviscid.  Numeri¬ 
cal  solution  of  the  Euler  equations,  though  expensive,  is  much  less  ex¬ 
pensive  than  solution  of  the  Navier-Stokes  equations  primarily  because 
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of  the  elimination  of  the  computational  work  and  computer  storage  asso¬ 
ciated  with  resolving  flow  variations  at  the  small  length  scale  at  which 
viscous  effects  are  important.  Time-accurate  numerical  solutions  of  the 
Euler  equations  about  three-dimensional  wings  in  oscillatory  motion  have 
been  previously  presented  by  Sankar,  et  al.,1  and  by  Salmond.2  Two-di¬ 
mensional  Euler  or  Navier-Stokes  calculations  about  oscillating  airfoils 
in  transonic  flow  have  been  presented  by  Magnus  and  Yoshihara,^  Steger 
and  Bailey,14  Chyu  and  Davis, ^  Janus, ^  and  Smith,  et  al.^  Time-accurate 
Euler  calculations  about  a  body  of  revolution  moving  away  from  a  flat 
plate  were  reported  in  Reference  8.  The  explicit  solution  algorithms 
used  in  References  3.  6,  7,  and  8  allowed  only  very  small  time  steps  to 
be  taken  because  of  stability  limitations.  Implicit  solution 
algorithms,  as  used  in  References  t,  2,  4 ,  and  5  typically  allow  much 
larger  time  steps  to  be  taken  while  maintaining  a  stable  solution,  but 
the  effect  of  such  large  time  steps  on  time-accuracy  is  unclear.  In 
Chapter  II  of  this  paper  the  Euler  equations  are  presented  and  trans¬ 
formed  from  rectangular  Cartesian  coordinates  into  general  time-depend¬ 
ent  curvilinear  coordinates  to  allow  for  dynamic  grids  that  can  follow 
the  motion  of  the  body.  In  Chapter  III  an  implicit  upwind  finite  volume 
scheme  is  presented  for  solving  these  equations,  and  several  variations 
of  the  scheme  are  applied  to  an  airfoil  oscillating  in  pitch  in  tran  - 
sonic  flow  to  study  the  effects  of  different  time  step  sizes,  different 
approximate  factorizations,  and  first-order  time-accuracy  versus 
second-order  time-accuracy  on  the  results. 

A  pacing  factor  for  the  ability  to  obtain  steady  or  unsteady 
numerical  solutions  of  the  Euler  equations  about  complex  geometries  is 
the  generation  of  a  mesh  of  points,  or  grid,  around  these  geometries. 


The  blocked  grid  approach  described  in  Chapter  IV  allows  much  greater 
flexibility  in  grid  generation  by  allowing  different  curvilinear  coordi¬ 
nate  transformations  to  be  used  for  different  parts  of  the  grid,  with 
only  the  requirement  that  all  grid  lines  be  continuous  from  one  block  to 
the  next  across  the  block-to-block  boundaries.  A  scheme  for  solving  the 
unsteady  Euler  equations  on  blocked  grids  by  advancing  the  solution  one 
time  step  in  each  block  sequentially  is  presented  in  Chapter  IV.  Tech¬ 
niques  for  transfer  of  information  across  block  boundaries  are  discussed 
along  with  analyses  of  error  associated  with  the  various  techniques. 
Steady  and  unsteady  calculations  on  blocked  grids  around  an  airfoil  are 
compared  with  calculations  on  equivalent  unblocked  grids  to  verify  the 
analyses  and  to  investigate  the  effect  of  block  boundaries  on  moving 
shock  waves . 

Three-dimensional  calculations  on  a  rectangular  planform  supercrit¬ 
ical  wing  of  aspect  ratio  two  are  presented  in  Chapter  V  for  both  steady 
flow  and  oscillatory  pitching  of  the  wing,  and  the  computed  results  are 
compared  with  the  experiment. 


CHAPTER  II 


THE  EULER  EQUATIONS 


The  three-dimensional,  time  dependent  Euler  equations  written  in 
3trong  conservation  law  form  are 
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Here  p  is  the  mass  density  of  the  compressible  fluid,  u,  0,  and  w 
are  the  velocity  components  in  the  x,y,  and  z  Cartesian  coordinate  di¬ 
rections,  respectively,  p  is  the  pressure,  and  e  the  total  specific  en¬ 
ergy  of  the  fluid.  The  relationship  between  e  and  p  is  the  result  of 
making  the  perfect  gas  assumption,  and  Y  is  the  ratio  of  specific  heats. 
Carets  on  the  quantities  in  Equation  (2-1)  are  used  to  indicate  that 
dimensional  quantities  are  being  used.  One  way  to  non-dimensionalize 
Equation  (2-1)  is  to  define  the  following  non-dimensional  quantities: 


where  tr  is  any  convenient  reference  length,  (^p./p*) 1  ,  the  speed 

of  sound  in  undisturbed  gas,  and  is  the  density  of  undisturbed  gas. 
After  some  manipulation,  the  non-dimensional  Euler  equations  may  be 
written  in  exactly  the  same  form  as  Equation  (2-1)  with  all  carets 
dropped. 

To  simplify  numerical  treatment  of  boundary  conditions  around  gen¬ 
eral  geometries  the  Euler  equations  are  recast  in  terms  of  general  boun¬ 
dary  conforming  curvilinear  coordinates.  The  curvilinear  coordinates 
are  defined  a3 

£  -  £(x,y,z,t) 

n  -  n(x,y,z,t)  (2-3) 

C  -  e(x,y,z,t) 
t  -  t  ( t ) 

Note  that  the  curvilinear  coordinates  are  time  dependent  thus  allowing 
constant  coordinate  surfaces  to  follow  the  unsteady  motion  of  the  bound¬ 
ary.  Also  note  that  the  time  coordinate  is  not  dependent  on  spatial 


location. 


Applying  the  transformation  (2-3)  to  the  non-dimensional  Euler 


equations  (Equation  (2-1)  with  all  carets  dropped)  gives 


where 


3«  3F  3C  4  3H  . 

TrTTTnlc"0 


(2-4) 


Q  -  J  pv 


puU+£xP 


F  -  tTJ  pvU*tyP 

pwU+Ezp 

U(e+p)-£tp 


puV+nxp 

tTJ  PVV*DyP 

pwV*nzp 

V(e*p)-ntp 


puW+;„p 


and  U-Exu*Eyv*E2w*et, 


and  V-nxu+nyv+nzw+nt , 


tTJ  pvW+CyP 
pwW+czp 
M(e*p)-C»p 


and  W»txu+ty v+^zw+tt . 


In  the  above  equation,  J  is  the  Jacobian  of  the  inverse  transformation, 
i.e.  3(x,y,z)/3U,n,c) .  and  is  given  by 


J  -  x^(ynzc  -  zny^)  "  yc(xnzc  -  znx?)  +  "  ynxc^' 

The  metric  quantities  are  given  by 

5X  "  J-1(ynz?  “  zny;^  nx  “  J_1(z5yc  “  y£izC^ 

Cy  »  J-1(znx^  -  xnzc)  ny  -  -  z^Xj.) 

£ z  »  J_1(xny[.  -  ynx;)  n2  -  J-1(x^y^  -  y^x^) 

Ct  -  (-xTSx  -  yT^y  -  zT£;2)Tt  nt  -  (~xTnx  -  yTny  -  zTnz)xt 

Cx  “  J_1(y5zn  “  z5yn) 

Cy  -  J-1(xnZj.  -  znx^) 

<z  “  J_1(Vn  "  yCxn} 

Ct  “  ^_xt^x  ~  yT?y  _  zx^z^Tt 


The  details  of  this  transformation  are  given  In  Appendix  A. 


CHAPTER  III 

IMPLICIT  SOLUTION  ALGORITHM 


3.1  FLUX  SPLITTING 

Hyperbolic  partial  differential  equations,  such  as  the  Euler  equa¬ 
tions,  are  characterized  by  the  existence  of  a  limited  domain  of  depend¬ 
ence.  The  solution  at  a  point  does  not  depend  on  every  other  point  in 
the  field;  this  means  that  information  travels  only  in  certain  charac¬ 
teristic  directions.  Numerical  schemes  intended  to  solve  hyperbolic 
equations  are  usually  enhanced  by  insuring  that  the  numerical  method 
propagates  Information  in  the  direction  specified  by  the  partial  differ¬ 
ential  equation.  This  can  be  done  by  using  an  upwind  method,  or  one  in 
which  the  difference  operator  is  taken  in  the  direction  from  which  the 


m 


V/, 

#•; 
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information  should  come.  Stability  properties  are  often  improved  by 
upwinding,  and  it  is  usually  unnecessary  to  add  smoothing  terms  or  arti¬ 
ficial  viscosity  to  an  upwind  method. 

The  three-dimensional  Euler  equations.  Equation  (2-4),  are  a  hypei 
bolic  system  of  five  equations  and  hence  have  five  characteristic  ve¬ 
locities  in  each  of  the  three  spatial  directions.  These  characteristic 
velocities  are  determined  from  the  quasillnear  form  of  Equation  (2-4), 
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where  the  matrices  A,  3,  and  C  are  given  by 
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The  eigenvalues  of  A  are  the  characteristic  velocities  in  the  £  direc¬ 
tion;  the  eigenvalues  of  B  are  the  characteristic  velocities  in  the  n 
direction;  and  similarly  the  eigenvalues  of  C  are  the  characteristic 
velocities  in  the  c  direction. 

Since  F,  G,  and  H  are  identical  except  that  where  £  appears  in  F,  n 
will  appear  in  G  and  £  will  appear  in  H,  extra  writing  can  be  avoided  by 
letting  K  represent  either  F,  G,  or  H,  when  k  represent  either  £,  n,  or 
C,  respectively.  Then  define 


(3-3) 


which  corresponds  to  A,  B,  or  C  depending  on  the  meaning  of  K.  The  gen¬ 
eral  flux  Jacobian  matrix  K  and  its  eigenvalues  are  derived  in  Reference 


6.  The  eigenvalues  of  the  matrix  K  are 


Ak  "  *  ^kxu  +  kyv  +  kzv  *  ^t^T  “  8ktx 


-  (Bk  +  c !  Vk  I  )tT  -  (8k  ♦  c(k^  ♦  ky  +  kl)W2)tT 


Ak  -  (8k  -  c I Vk I ) tT 
where  c  is  the  speed  of  sound. 

It  is  possible  to  split  tne  flux  vector  K  into  three  parts,  one 
corresponding  to  each  of  the  distinct  eigenvalues  of  K  given  above.  (For 
the  details  of  this  splitting,  see  Reference  6.)  The  flux  vector  K  is 
then  written  as 


K  ’  XkK1  +  AkKi4  *  Me 
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The  sign  of  in  Equation  (3-10  determines  from  which  direction  infor 
mation  should  be  used  to  determine  the  corresponding  portion  of  flux 


K,  ,  for  2.-1  ,  4,  and  5. 
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The  discussion  above  provides  a  rationale  for  writing  the  flux  vec- 

L  ft 

tor  K  as  the  sum  of  two  vectors,  K  and  K  , 

K  -  KL  ♦  KR.  (3-7) 

f  O 

KL'  is  associated  with  the  eigenvalues  that  have  positive  signs,  and  K 

is  associated  with  the  eigenvalues  that  have  negative  signs.  The  spe- 

I  H 

cific  method  used  to  calculate  K  and  K  will  be  discussed  in  the  next 
section. 


3.2  Discretization 

A  finite  volume  discretization  of  Equation  (3~1)  balances  the  in¬ 
crease  of  the  conserved  quantity  in  a  computational  cell,  or  volume, 
with  the  flux  of  the  quantity  through  the  surface  of  the  cell.  Figure 
(3“1 )  depicts  a  portion  of  the  computational  domain  with  a  typical  cell 
labeled.  Assuming  the  dependent  variables  are  constant  in  the  interior 
of  cell  (i,  j,  k),  and  that  the  flux  vectors  F,  G,  and  H  are  constant 
over  the  constant  constant  n,  and  constant  c  surfaces  of  the  cell, 
respectively,  an  implicit  discretization  of  Equation  (2-*0  is  (Reference 
9) 

(QnM  ,  -  Qn  .  ,  )AU nAc  *  (Fn\  -  Fn*;  )UndClT  .3-3) 

,J’*  i--,j,k 
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This  can  be  written  as 


AtU^F0*1  ♦  5nGn4-1  *  6cHn+1)  -  0 


(3-9) 


where  <5^,  for  example,  is  defined  by 


VCi,J.K  “  IT-(f  i  "  F  _l  K 
i * 2 • J  ^"2 » j 


and  6n  and  6,  are  defined  analogously.  Letting  F^,F,:i,Gu,G“,H“,  and  H" 
be  the  split  flux  vectors  for  F,  G,  and  H  as  given  by  Equation  (3_7),  an 
implicit  split  flux  discretization  is  given  by 

AQn  +  At[6^(F“  *  Fa)n+1  *  5n(GL  *  G*)"*1  *  6 j- ( H“  ♦  ]  -  3  (3*10) 

The  fluxes  are  nonlinear  functions  of  the  dependent  variables,  and 
must  be  linearized  to  obtain  a  system  of  equations  that  can  easily  be 
solved  for  AQn.  Beam  and  Warming10  and  3riley  and  McDonald11  have  used 


the  linearization 


Kn>1  -  Kn  ♦  RnAQn  ♦  0(At2), 


(3-n) 


where  -  (3K/3Q)n  for  a  typical  flux  terra  K.  The  flux  Jacobian  matri¬ 
ces  required  to  linearize  Equation  (3“)0)  are 
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and  are  derived  in  Appendix  B.  Using  this  linearization,  Equation 
( 3~ 10)  becomes 


I  ♦  At  ( A*1  ♦  5*  A*  +  6 1B^‘  ♦  6  S?  ♦  6iC*  ♦  61C?)]AQn 
5  £  n  n  c  c 


*At  ( <S ®Fn  ♦  6eGn  +  6®Hn) , 
E  n  5 


(3-13) 


where  a  distinction  has  been  made  between  the  implicit  spatial  differ¬ 
ence  operators  and  the  explicit  spatial  difference  operators  by  using 
superscripts  i  and  e  respectively.  The  dots  are  used  to  indicate  that 


the  difference  operators  apply  to  the  product  of  the  Jacobian  matrices 
with  AQn.  If  6®,  6®,  and  6®  are  chosen  such  that 
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and  6^ ,  5^,  and  <5 *  are  chosen  such  that 
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then  substitution  into  Equation  ( 3— 1 3 )  yields 
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*  0(At  ) )  -  -[|£  ♦  |5  ♦  ♦  0(A£2,An2,AC2)], 


after  dividing  through  by  At  and  using  the  fact  that 


AQn  -  Qn+1  -  Qn  -  At  |S  ♦  0(At2). 


(3-17) 


Using  Equation  (3-11)  in  (3-16)  results  in 
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Therefore  Equations  (3-1*0  and  { 3— 15)  give  a  method  that  is  second-order 
accurate  in  space  if  A£ ,  An,  and  Ac  are  O(At)  and  first-order  accurate 
in  time  when  used  with  Equation  (3-13).  There  are  at  least  two  simple 
variations  of  this  algorithm  that  will  produce  second-order  accuracy  in 


Before  proceeding  to  a  discussion  of  these  second-order  time-accu¬ 
rate  methods,  simplify  the  problem  by  considering  only  one  spatial  di¬ 
mension.  Also,  observe  from  the  analysis  above  that  carrying  the  split 
flux  vectors  had  no  effect  on  the  accuracy  analysis.  The  fact  that  this 
is  a  split-flux  method  only  enters  the  problem  when  actually  defining 
the  difference  operators  to  satisfy  the  conditions  given  by  Equations 
(3-1*0  and  (3“15),  which  has  not  yet  been  done. 

The  simplified  problem  that  will  now  be  considered  is 
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3ecause  the  time  accuracy  of  the  method  is  independent  of  how  the  split 
fluxes  are  spatially  differenced,  any  standard  time  discretization  can 
be  applied.  A  broad  class  of  schemes  is  given  by  Beam  and  Warming's 
formula1 2 


U  ♦  4Qn"'  *  one  -  J,  - 


(3-20) 


The  scheme  given  by  (3-13)  is  the  backward  Euler  differencing  scheme 
obtained  from  (3-20)  by  setting  9-1  and  v0.  By  taking  9-1/2,  '>-0,  a 
trapezoidal  scheme  is  obtained  that  is  second  order  accurate.  A  three 
point  backward  scheme  is  obtained  by  setting  9-1,  H/-1/2,  giving 
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Spatially  differencing  Equation  (3-21)  with 
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0(AtAx2, At2Ax,At3) .  (3“23) 


For  the  3ake  of  completeness,  the  trapezoidal  scheme  for  Equation  (3—19) 
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and  the  backward  Euler  scheme  written  for  Equation  (3—19)  is 


[i  *  At6iA^]AQn  -  -At6eF  *  0(AtAx2,At2Ax, At2) . 
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Although  the  truncation  error  has  been  determined  without  fully 
specifying  the  spatial  difference  operators,  other  properties  of  the 

i  p 

various  methods  cannot  be  determined  until  5  and  5  are  completely  de¬ 
fined.  To  be  consistent  with  the  objective  of  deriving  a  finite  volume 
method,  61  and  5e  should  be  central  differences  of  quantities  at  oppos¬ 
ing  faces  of  a  computational  cell.  Dependent  variable  values  are  only 
stored  at  cell  centers,  so  some  method  must  be  specified  for  determining 
values  at  cell  faces.  Also,  the  difference  operators  must  respect  the 

I  D 

direction  of  information  transfer  and  difference  F  and  F  in  their  re¬ 
spective  upwind  directions.  These  two  requirements  are  satisfied  by 
defining 


6epj  "  Ix[(F"L  t  ~  )  -  (FR  -  FR  )]  (3-26) 

where  Fj>1/2  depends  only  on  Information  to  the  left  of  point  j+1/2, 

I  R 

i.e.  Fj+i/2  depends  on  Qj,Qj_1t...,  and  F 1/2  dePends  only  on  informa¬ 
tion  to  the  right  of  point  J-H/2,  i.e.  Fj  +  1/,2  dePends  on  Qj  +  i>  Oj^***- 
One  method  to  satisfy  the  accuracy  requirement  given  by  Equation 
(3-22)  for  6e  is  to  evaluate  FL  and  FR  by  two  point  flux  extrapolations 
given  by 
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If  the  eigenvalues  do  not  change  sign  from  the  j -1/2  face  to  the  j  +  1/2 
face,  substitution  of  Equations  (3"27)  into  Equation  (3-26)  gives  either 
a  second  order  backward  or  a  second  order  forward  difference  approxima¬ 
tion  to  each  split  flux  component  yielding 
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(3-28) 


This  is  true  only  if  the  eigenvalues  are  of  the  same  sign  on  both  faces 

of  the  cell.  Also  notice  that  there  is  some  ambiguity  in  using  Equation 

0 

(3~27)  because  we  use  Xj+1/2  which  i3  at  the  cel1  face,  but  have  not 
specified  how  to  obtain  these  quantities. 

Another  alternative  for  evaluating  F  and  Fn  30  that  Equation 
(3-26)  is  a  second  order  accurate  first  derivative  of  F  is  to  evaluate 
FL  using  a  two  point  dependent  variable  extrapolation  from  the  left  and 

p 

F  using  a  two  point  dependent  variable  extrapolation  from  the  right. 
This  i3  the  method  of  Reference  8  and  is  used  for  all  calculations  re¬ 
ported  here.  For  a  particular  cell  face,  the  scheme  reported  here  uses 
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to  calculate  a  3et  of  left  eigenvalues,  X*(QL),  and  left  split  fluxes, 
F1(Qt"'»,  and  uses 
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to  calculate  a  set  of  right  eigenvalues,  X  (Q  ),  and  right  split  fluxes, 


F4(Qr).  The  flux  at  the  face  is  then  set  to 

FL  -  l  i  (X*(QL)  ♦  lxt(QL)l)F,(QL)  (3-31a) 

l  2 

and 

FR  -  l  i  (Xl(QR)  -  lx1(QR)l)Fi(QR)  (3-31 b ) 

l 


Note  that  if  a  left  and  right  eigenvalue  have  different  signs,  then  the 
corresponding  split  flux  may  be  either  summed  from  both  sides,  if  the 
3ign  changes  from  positive  to  negative,  or  from  neither  side,  if  the 
sign  changes  from  negative  to  positive. 

To  determine  the  truncation  error  of  this  scheme  consider  the  sca¬ 
lar  version  of  Equation  ( 3“ 1 9 )  and  assume  3F/3Q>0  for  all  x.  Then  FR-0 
and  Equation  (3-26)  would  become 

*@FJ  "  Sx^^J  "  “  7  QJ-2^  (3-32) 


Letting  AQj-Qj-Qj_,  and  using  Taylor  series  expansion  about  Qj 

F(|  Qj  -  \  Qj.,)  -  F(Qj  *  j  »«j)  -  F!0j>  *  j  «j  lo'j 
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Substituting  Equations  (3-33)  and  (3-3*0  in  Equation  (3~32),  followed  by 
expanding  all  dependent  variable  values  about  Xj  and  then  simplifying 


yields 
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3y  comparing  Equations  (3-28)  and  ( 3— 35 )  it  is  clear  that  extrapolating 
dependent  variables  is  different  from  extrapolating  flux  values  only  in 


the  higher  order  terms. 

The  implicit  difference  operation  61  need  only  be  first  order  accu¬ 


rate  and  hence  it  is  evaluated  using  one  point  dependent  variable  ex¬ 
trapolation  giving 


6 1  (AAQ)j  -  (aL(qJ)AQJ  -  AL(Qj_1)dQ"_1) 


♦  ( AR(Qj  +  1  )4Qj  +  1  -  AR(Qjn)AQj)  ^ 


(3-36) 


Metric  quantities  involved  in  the  evaluation  of  the  Jacobian  matrices 
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are  taken  from  the  appropriate  face,  so  that  Al(Qj)  uses  metrics  from 
j+1/2  and  AR(Qj)  uses  metrics  from  j-1/2  in  Equation  (3-36). 

For  dynamic  grid  calculations,  metric  quantities  including  the  cell 
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volumes  change  during  the  solution  process.  For  this  study,  all  spatial 
metrics  are  evaluated  at  time  level  n  to  be  used  in  the  calculation  of 
AQn.  The  time  metrics  at  time  level  n  are  evaluated  from  the  change  in 
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3.3  Factorizations 


The  finite  volume  scheme  described  by  Equation  (3-13)  is  not  prac¬ 
tical  to  use  because  the  system  of  algebraic  equations  that  are  gener¬ 
ated  have  a  very  large  bandwidth.  There  are  many  possible  approximate 
factorizations  of  this  method  that  can  be  used  to  make  the  solution  pro¬ 
cess  easier,  as  noted  by  Steger  and  Warming,111  for  example.  Several 

Q 

factorizations  have  been  recently  applied  by  Whitfield,  and  by  Ander¬ 
son,  Thomas,  and  Whitfield. 1  One  scheme  that  has  been  used  is  the 
six-factor  scheme  given  by 

[l+AT6^A-][r+AT6^A?][l+AT6^B*][l+AT<5^8*][l+AT5^C-][l*AT6^C-]AQn  (3-39) 

-  -AxRn, 

where 

Rn-6|Fn  *  «« Gn  *  5«Hn 

The  solution  of  this  scheme  consists  of  successive  solution  of  six  block 
bidiagonal  systems  of  equations  which  are 


L 


AtA^A- ]x’ 

-  -AiRn 

(3-^0a) 

R 

AtA|A-]x2 

-  X1 

(3-^0b) 

L 

AtA^B* ]x3 

■  X2 

(3-^0c) 

R 

At A^B* ]XU 

-  X3 

(3-U0d) 

,v!v 


(3-^Oe) 


[i  ♦  ATijc^jx5  -  X4 
R 

[i  ♦  At6^C • ] AQn  -  X5  (3-^Of) 

The  bidiagonal  solution  can  be  vectorized  in  planes  perpendicular  to 
lines  along  which  backward  or  forward  substitution  is  taking  place.  For 
example,  the  solution  of  (3“*<0a),  which  is  lower  bidiagonal,  is  obtained 
by  forward  substituting  along  £  lines.  In  this  case  the  vectori zation 
direction  is  the  C-constant  plane.  Another  factorization  is 

[l+AT^A-*5^A?)][l*AT(5^B**6^B-)][l*AT(6^C*+5£c*)]AQn  -  -AtRn.  (3"^) 

The  solution  is  obtained  to  this  three-factor  scheme  by  successive  solu¬ 
tion  of  three  block  tridiagonal  systems  of  equations,  and  can  be  vec¬ 
torized  similarly  to  the  six-factor  scheme.  This  method  requires  half 
the  solution  passes  of  the  six-factor  scheme,  but  each  pass  is  more 
costly  since  it  requires  solution  of  a  block  tridiagonal  system.  This 
method  also  requires  more  storage  than  the  six-factor  because  two  sets 
of  flux  Jacobians  are  required  in  memory  at  the  same  time,  while  the 
3ix-factor  only  needs  one  at  a  time.  Yet  another  factorization  is  given 
by 

[  I  ♦At  ( 5^A^+6^B*'  +6 ] [ I+At  ( 5^  A  - *i5^B* +  )  ]  AQn  -  -AiRn  (3~^2) 

This  two-factor  scheme  is  very  attractive  in  that  it  consists  of  solu¬ 
tion  of  a  sparse  block  lower  triangular  system  followed  by  solution  of  a 
sparse  block  upper  triangular  system  of  equations  given  by 


22 


*  •  k. 


[i  +  Ax  (5^A*  ♦  <5^B*  ♦  &lQC‘]X]  -  -AxRn 
[i  ♦  At (6pA^  ♦  6*B*  +  5^C*]AQn  -  X1 . 


(3-^ 3a) 


( 3-^3t> ) 


Solution  of  (3_i*3a)  is  done  by  a  simple  forward  substitution  and  solu¬ 
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tion  of  (3-1*3b)  by  a  simple  back  substitution.  As  noted  in  Reference 
15,  the  forward  and  backward  substitution  steps  cannot  be  vectorized 


along  any  coordinate  line  or  constant  coordinate  plane.  This  is  due  to 
the  fact  that  xj  j  k  in  ( 3— ^ 3a )  depends  on  x|_1  j  k,  X ^ » j  —  i  i<  and 
xi  j  k-1’  ‘However,  Janus1**  has  subsequently  derived  an  algorithm  for 
vectorizing  the  substitution  step  by  operating  on  all  points  in  a  diago¬ 


nal  plane  in  computational  space  simultaneously.  A  diagonal  plane  is 


one  on  which  i+j+k  ■  constant.  The  solution  of  Equation  ( 3“ ^ 3 )  is  a 
pass  through  computational  space  starting  at  the  lower  corner,  i.e. 


y.-i 
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( i  ,j  ,k)-( 1 , 1  , 1 ) ,  and  sweeping  through  by  diagonal  planes  until  the  upper 
corner,  ( i ,J ,k )-( imax,  Jmax.  kmax)  ls  reached-  The  second  pass,  given 
by  Equation  ( 3~**3t> )  starts  at  the  upper  corner  and  sweeps  back  to  the 
lower  corner.  Using- this  technique  the  solution  of  Equation  (3_1*3)  is 
fully  vectorizable .  One  disadvantage  of  this  two-factor  algorithm  is 
that  three  sets  of  Jacobian  matrices  need  to  be  stored  simultaneously. 

A  stability  analysis  for  the  six-factor,  three-factor  and  two  fac¬ 
tor  schemes  was  presented  by  Anderson,  et  al.,  in  Reference  15.  The 
analysis  indicated  that  the  six-factor  scheme  had  the  poorest  stability 
properties  and  was  stable  up  to  a  Courant  number  of  about  10.  The 
two-factor  scheme  was  least  sensitive  to  Courant  number  and  was  stable 
for  all  Courant  numbers  up  to  35,  which  was  the  maximum  investigated. 
The  present  study  includes  calculations  with  both  the  six-factor  and  the 
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two-factor  methods  which  tend  to  verify  the  superiority  of  the  two  fac¬ 
tor  results  for  higher  Courant  numbers.  These  results  will  be  presented 
in  Section  3-5. 

3.4  Boundary  Conditions 

Characteristic  variable  boundary  conditions  are  derived  in  Refer¬ 
ence  17  for  stationary  grids  and  Reference  8  for  dynamic  grids.  A3  in 
these  references,  a  layer  of  points  outside  the  boundaries,  called  phan¬ 
tom  points,  are  used  here  to  implement  the  characteristic  variable 
boundary  conditions.  Because  there  is  only  one  point  outside  the  bound¬ 
ary,  the  computation  of  fluxes  at  boundary  faces  is  modified  to  use  one 
point  dependent  variable  extrapolation  at  boundary  faces  both  from  in¬ 
side  and  outside  the  domain.  The  change  in  dependent  variables,  AQn,  is 
taken  as  zero  at  phantom  points.  The  result  of  this  boundary  treatment 
is  that  the  boundary  conditions  are  only  first  order  accurate.  The  as¬ 
sumption  is  made  that  this  will  not  affect  global  accuracy  of  the  second 
order  scheme. 

3.5  Time  Accuracy 

In  this  section  the  time-accuracy  of  the  split  flux  algorithms  pre¬ 
viously  discussed  will  be  examined  by  comparing  computed  unsteady 
results  for  different  factorizations,  different  time  discretizations, 
and  different  time  3tep  sizes.  All  calculations  are  for  a  NACA0012  air¬ 
foil  oscillating  in  pitch  about  its  quarter  chord  point  at  a  Mach  number 
of  0.755  and  a  reduced  frequency  of  0.1628,  where  reduced  frequency,  or 
Strouhal  number,  is  defined  as 


with  w  -  oscillatory  frequency  of  airfoil  motion 
V,,  »  free-stream  velocity 
$  »  chord  length  of  airfoil 

The  mean  angle  of  attack,  aQ,  of  the  airfoil  was  0.016  degrees  and  the 
unsteady  angle  of  attack  amplitude,  a-, ,  was  2.51  degrees.  The  angle  of 
attack  history  of  the  airfoil  is  given  by 

a(t)  -  aQ  +  a1  sin(2it),  (3-iJ5) 

which  in  nondimens ional  variables  is  given  by 

a(t)  -  aQ  ♦  a1  sin(Ma>kt)  (3_u6) 

1  ft 

Experimental  data  is  available  for  this  case  from  Landon,  °  and  a  com¬ 
parison  of  computed  results  with  the  experimental  data  will  also  be  pre¬ 
sented. 

The  grid  used  for  these  calculations  was  a  'C'  grid  with  221  points 
on  the  airfoil  and  in  the  wake  and  20  points  from  the  airfoil  to  the 
far-field  boundaries.  The  far-field  boundaries  were  placed  25  chords 
upstream  and  30  chords  above,  below,  and  downstream  of  the  airfoil. 
There  were  1 6 1  points  actually  on  the  airfoil  surface.  The  grid  was  gen¬ 
erated  using  the  transfinite  interpolation  option  of  the  grid  code  writ¬ 
ten  by  Thompson.  Figures  (3-2a)  and  (3~2b)  show  the  entire  grid  and  a 
cl03eup  of  the  airfoil.  The  maximum  spacing  between  points  on  the  air¬ 
foil  surface  is  1.9%  of  the  chordlength  and  occurs  at  the  60%  chord  po¬ 
sition.  From  experience  with  steady  calculations  on  this  grid  it  has 
been  found  that  this  grid  resolution  is  insufficient  to  resolve  weak 


shocks,  but  clustering  the  grid  at  shock  locations  is  Impractical  for 
unsteady  calculations  in  which  the  shock  position  changes  by  20i  of  the 
chord  during  the  calculation.  Solution  adaptive  gridding  offers  a  solu¬ 
tion  to  this  problem,  but  is  beyond  the  scope  of  this  investigation. 

The  calculations  were  started  with  Initial  conditions  of  free 
stream  values  everywhere,  and  then  500  iterations  of  local  time  stepping 
with  CFL  (Equation  ( 3— 38 ) )  equal  to  15  were  used  to  establish  the  steady 
flow  field  before  the  unsteady  motion  was  started.  The  local  time  step¬ 
ping  was  replaced  by  a  fixed  time  step,  the  same  at  all  points,  when  the 
unsteady  calculation  began.  The  unsteady  calculations  were  continued 
until  the  airfoil  had  gone  through  four  cycles  of  pitching  motion  to 
allow  unsteady  transients  to  decay.  Airfoil  motion  for  this  and  all 
cases  presented  here  was  Implemented  by  rotating  the  entire  grid  as  a 
rigid  body  about  the  pitch  axis. 

Results  from  the  two-pass  algorithm  applied  to  the  oscillating 
NACA001 2  case  described  above  are  presented  in  Figures  (3~3)  through 
(3-5).  Figures  (3~3a  thru  h)  show  the  pressure  coefficients  on  the  air¬ 
foil  surface  at  eight  time  instants  during  the  fourth  cycle  of  oscillat¬ 
ing  pitch.  Three  time  step  sizes  are  shown  using  the 
second-order  time-accurate  three  point  backward  scheme  given  by  Equa¬ 
tion  (3-23),  and  two  time  step  sizes  using  the  first-order  time-accurate 
backward  Euler  scheme  given  by  Equation  (3”25).  In  Figure  (3~3a)  the 
airfoil  is  25°  into  the  cycle,  and  although  the  angle  of  attack  i3  in¬ 
creasing  through  1.09°,  there  is  a  remnant  of  the  lower  surface  shock 
wave  moving  forward  and  weakening.  This  phenomenon,  and  the  similar  oc¬ 
currence  involving  the  collapse  of  the  upper  surface  shock  when  the  an- 


gle  of  is  decreasing  (Figure  3~3e),  is  where  the  variation 
between  the  different  methods  and  time  steps  is  most  visible.  The  ex¬ 
tremes  in  the  shock  locations  given  by  the  various  curves  are  always 
given  by  the  0.01  time  step  case  and  the  0.20  time  step  case.  This  is 
true  at  each  instant  plotted  in  Figures  (3-3a-h).  In  all  cases,  the 
shock  wave  given  by  the  larger  time  step  lags  the  motion  of  the  shock 
wave  given  by  the  smaller  time  step.  Taking  the  shock  location  as  the 
position  where  the  coefficient  of  pressure  curve  crosses  the  critical 
pressure  coefficient  value,  Figure  (3-4)  shows  the  shock  location  for 
all  the  cases  shown  in  Figure  (3_3a).  Notice  from  Figure  (3-4)  that 
decreasing  the  time  step  an  order  of  magnitude  from  At  -  0.10  to  At  - 
0.01  changed  the  shock  location  by  approximately  the  same  amount  as  re¬ 
ducing  the  time  step  from  At  -  0.20  to  0.10.  This  indicates  that  reduc¬ 
tion  of  the  time  step  below  At »0.01  should  not  be  expected  to  change  the 
shock  location  appreciably.  Returning  to  Figure  (3~3),  another  differ¬ 
ence  between  results  using  different  time  step  sizes  is  apparent.  For 
those  cases  with  relatively  large  shock  wave  speeds,  the  larger  time 
steps  show  more  compression  downstream  of  the  shock  and  a  large  wave¬ 
length  'ringing'  behavior,  but  in  parts  (c)  and  (g),  where  the  shock  has 
nearly  reached  the  point  of  maximum  downstream  travel,  the  smaller  time 
step  produces  a  higher  pressure  coefficient  behind  the  shock,  with  no 
ringing,  however.  The  At-0.01  case  also  generally  produces  slightly 
sharper  shocks  than  the  larger  time  steps. 

Comparing  the  first-order  time-accurate  results  with  the  second- 
order  time-accurate  results  in  Figure  (3”3)  shows  little  difference.  The 
first-order  At-0.05  pressures  are  between  the  second-order  At-0.01  and 


'y  j _--r»  I1- '.'/r ’ts,  '.'v '.'■  *'' 1 


A 


second-order  At-0.10  pressures,  while  the  first-order  At-0.10  results 
are  slightly  more  different  from  the  At-0.01  results  than  are  the  sec¬ 
ond-order  At-0.10  results.  Since  the  second-order  At-0.01  case  can  be 
assumed  closer  to  the  correct  values,  this  indicates  that  the  second- 
order  method  gives  slightly  better  results  for  large  time  step  sizes, 
but  no  appreciable  difference  is  shown  for  this  test  case. 

The  difference  in  pressure  coefficients  for  the  different  time  step 
sizes  and  methods  discussed  above  make  relatively  little  difference  to 
the  lift  coefficient  history  of  the  airfoil,  as  shown  in  Figure  (3~5a), 
but  make  more  difference  to  the  moment  coefficient  history,  as  shown  in 
Figure  (3~5b).  The  ragged  peaks  are  due  to  not  calculating  lift  and 
moment  at  the  same  time  in  each  cycle.  One  reason  for  the  small  differ¬ 
ence  in  lift  is  that  the  increased  area  under  the  pressure  coefficient 
curve  due  to  the  lagging  shock  position  of  the  larger  time  step  cases  is 
counterbalanced  by  the  over compress ion  behind  the  shock.  This  behavior 
imposes  an  additional  moment  on  the  airfoil,  however,  and  is  responsible 
for  the  larger  amplitudes  of  the  At-0.10  and  At-0.20  moment  coefficient 
histories . 

The  discussion  above  has  dealt  with  the  time  step  sizes  used  for 
different  calculations.  Time  step  size  by  itself  has  little  meaning  for 
application  to  general  problems  that  may  have  different  arids  and  flow 
conditions.  A  parameter  that  is  much  more  meaningful  for  numerical  so¬ 
lution  of  hyperbolic  systems  is  the  Courant  number.  The  Courant  number 
is  the  number  of  computational  cells  that  a  wave  will  propagate  through 
in  one  time  step.  For  a  given  time  3tep,  At,  Equations  (3-38a,b)  give 
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a  value  for  the  CFL  parameter  that  will  be  used  as  a  measure  of  Ccurant 


number  for  the  three-dimensional  problem. 


CFL 


Ax 


max 


max  lx*' I  max  I  A*'! 
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(3-^7) 


The  maximum  CFL  for  the  Ax-0.01  case  was  approximately  10,  the  maximum 
CFL  for  the  Ax-0.10  case  was  about  100,  and  for  the  Ax-0.20  case  the 
maximum  CFL  was  about  200.  The  maximum  CFL  occurs  in  the  first  cell  off 
the  airfoil  at  about  2.5X  chords  downstream  of  the  leading  edge.  At  the 
25S  chord  location,  the  CFL  is  about  half  its  maximum  value,  and  at  the 
trailing  edge  it  is  one-fourth  its  maximum.  The  local  CFL  decreases 
quickly  as  the  cells  become  larger  away  from  the  airfoil.  The  maximum 
CFL  in  the  second  row  of  points  off  the  airfoil  surface  is  only  75X  of 
the  maximum  CFL  in  the  first  row. 

All  results  presented  thus  far  were  obtained  using  the  two-pass,  or 
two-factor,  method.  The  six-pa33,  or  six-factor,  method  has  also  been 
applied  to  the  oscillating  NACA0012  airfoil.  Results  obtained  for 
Ax-0.01  and  Ax-0.05  are  similar  to  those  obtained  from  the  two-pass 
method,  but  for  the  Ax-0. 10  case  the  pressure  distribution  displayed 
oscillations  as  illustrated  in  Figure  (3~6a,b).  As  previously  dis¬ 
cussed,  the  maximum  CFL  for  this  case  is  near  100  and  occurs  not  at  the 
leading  edge  but  slightly  downstream  of  it  at  approximately  the  location 
where  those  oscillations  occur.  The  oscillations  are  present  on  both 
the  upper  and  lower  surface,  but  are  much  larger  in  magnitude  on  the 
lower  surface  throughout  the  cyclic  motion  of  the  airfoil,  even  though 
the  local  CFL  on  the  upper  surface  at  times  exceeds  that  of  the  lower 


surface  as  In  Figure  (3“6a).  The  reason  for  this  is  not  understood,  but 
may  be  due  to  the  asymmetry  of  the  six-pass  method.  These  results  verify 
the  trend  predicted  by  the  stability  analysis  of  Anderson,  et  al.1-’ 
discussed  earlier.  It  is  apparent  from  the  figure  that  in  regions  away 
from  the  oscillations,  the  six-pass  results  match  the  two-pass  results 
very  well. 

The  trapezoidal  scheme  given  by  Equation  ( 3— 24 )  was  attempted  for 
this  problem,  but  it  was  found  to  be  unstable  for  all  tested  values  of 
CFL  greater  than  one  using  local  time  stepping,  and  no  full  calculations 
were  performed. 

A  comparison  of  computed  pressure  coefficients  for  the  pitching 
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NACA001 2  with  the  experimental  results  of  Landon  is  given  in  Figure 
(3-7).  Figures  (3"8a,b)  are  plots  of  the  experimental  versus  the  com¬ 
puted  lift  and  moment  coefficient  histories.  Computed  results  are  plot¬ 
ted  for  both  the  nominal  mean  angle  of  attack,  Oq-0.016°,  and  for  a  mean 
angle  of  attack  of  a0-0.375°.  The  use  of  a  mean  angle  of  attack  dif¬ 
ferent  from  the  near  zero  angle  of  attack  reported  in  Reference  18  is 
motivated  by  the  differences  in  the  experimental  Cp  curves  at  points 
that  are  180°  apart  in  the  oscillatory  motion  of  the  airfoil  (compare 
Figures  (3~7b)  and  ( 3_ 7 f ) ) ,  and  by  the  obvious  bias  towards  positive 
lift  coefficients  shown  in  Figure  (3~8a).  The  specific  value  used  for 
aQ  was  obtained  by  using  the  slope  of  the  unsteady  lift  versus  angle  of 
attack  curve  to  determine  the  change  in  c»q  required  to  match  the  experi¬ 
mental  peak  lift  values.  The  a0-0.375°  calculation  matches  the  experi¬ 
ment  much  more  closely  than  the  nominal.  In  particular,  the 
experimental  data  in  Figure  (3"7a)  shows  an  already  developed  upper  sur- 


face  shock  in  addition  to  the  decaying  lower  surface  shock.  The 
ag-0.0l6°  case  gives  the  lower  surface  shock  too  far  downstream  and  too 
strong,  while  the  upper  surface  shock  has  not  yet  been  resolved.  The 
remaining  parts  of  Figure  (3-7)  also  illustrate  better  agreement  for  the 
a0«0.375°  case.  The  a0-0.375°  case  matches  the  experimental  lift  his¬ 
tory  better  than  the  ag-O.160  case  in  Figure  (3“8a),  but  both  computa¬ 
tions  fail  to  match  the  negative  portion  of  the  moment  coefficient 


history  in  Figure  ( 3— 8b ) . 
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CHAPTER  IV 


EULER  SOLUTIONS  ON  BLOCKED  GRIDS 


M.l  Basic  Concepts 


Various  approaches  have  been  presented  for  making  the  task  of  grid 
generation  about  complex  geometries  more  simple,  such  as  the  Chimera 
grid  scheme  of  Dougherty,  et  al.,1^  and  Benek,  et  al.20  Figure(4-ia) 
shows  a  schematic  of  a  Chimera  grid.  The  Chimera  approach  to  gridding 
provides  for  multiple  bodies  that  are  gridded  independently  and  communi¬ 
cate  between  grids  by  interpolating  in  regions  of  overlap.  The  patched 
grid  approach  of  Rai  requires  that  grids  meet  along  a  common  boundary 
as  shown  in  Figure  ( 4— 1b),  but  allows  for  discontinuity  in  the  grid 
across  the  boundary.  The  blocked  grid  scheme  advocated  by  Thompson, 
Figure  (*t-lc),  is  a  subset  of  patched  grids  that  assumes  complete  conti¬ 
nuity  (or  at  least  slope  continuity)  of  the  grid  across  block 
boundaries.  This  scheme  makes  the  task  of  information  transfer  across 
block  boundaries  much  simpler  than  when  grid  lines  are  not  continued 
into  adjacent  blocks  because  complex  interpolation  is  not  required  at 
grid  boundaries.  Though  blocked  grids  require  matched  Interface  grid 
boundaries  that  are  not  required  by  patched  grids,  blocked  grids  still 
allow  very  complicated  three-dimensional  configurations  to  be  easily 
gridded.  The  grid  used  by  Lasinski,  et  al.2^  to  calculate  flow  about  a 
tri-element  "augmentor  wing"  airfoil  is  an  example  of  a  blocked  grid 
used  to  model  a  complex  three-body  problem  in  two  dimensions.  Hessenius 
and  Pulliam21*  presented  steady  inviscid  two-dimensional  calculations  and 
unsteady  one-dinensional  calculations  on  blocked  grids.  An  implicit 
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central  difference  method  was  used  at  field  points  and  flux  splitting 

was  used  at  the  block-to-block  boundaries.  The  blocked  grid  concept  has 
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been  applied  to  three-dimensional  problems  by  Weatherill  and  Forsey 

and  by  Miki  and  Takagi.  The  common  theme  of  these  gridding  methods  is 

the  use  of  different  grids,  or  coordinates  systems,  in  different 

regions,  thereby  allowing  the  most  appropriate  grid  structure  to  be  used 

in  each  region.  The  primary  conceptual  difference  between  the  Chimera 

grids,  patched  grids,  and  blocked  grids  is  the  degree  of  independence 

between  the  various  subgrids.  The  idea  of  using  different  coordinates 
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systems  in  different  regions  is  not  new.  In  1969,  Thoman  and  Szewczyk 
calculated  viscous  flow  over  a  circular  cylinder  by  embedding  a  cylin¬ 
drical  coordinate  system  inside  a  rectangular  mesh.  The  cylindrical 
coordinates  allowed  better  resolution  near  the  cylinder,  while  the  rec¬ 
tangular  mesh  was  more  appropriate  for  the  wake  region. 

The  Euler  code  presented  here  will  accept  any  arrangement  of  arbi¬ 
trarily  sized  blocks.  Only  one  block  at  a  time  is  stored  in  memory 
while  all  other  blocks  are  stored  in  high  speed  secondary  memory  that  is 
accessed  like  a  disk.  Dynamic  memory  management  is  used  to  adjust  the 
amount  of  memory  requested  while  working  on  a  particular  block  to  equal 
the  amount  of  memory  needed.  Details  of  the  block  handling  and 
block-to-block  information  transfer  techniques  and  an  analysis  of  the 
error  introduced  at  block  boundaries  ill  also  be  presented  in  this 
chapter.  Calculations  for  transonic  flow  over  an  airfoil  using  blocked 
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II.  Note  that  I  ♦  AtLt  j-  and  I  +  AtLjj  xi  are  lower  block  triangular 
matrices  and  that  I  +  AtUt  j  and  I  +  AtUjj  n  are  uPPer  block  triangular 
matrices.  The  rectangular  matrix  AtLjj.j  contains  coefficients  of  quan¬ 
tities  in  block  I  that  are  required  to  make  spatial  difference  expres¬ 
sions  for  cells  in  block  II  complete.  Similarly,  contains 
coefficients  of  quantities  in  block  II  that  are  required  to  make  differ¬ 
ence  expressions  for  cells  in  block  I  complete.  Although  Equation 
(3-^3)  uses  backward  Euler  time  differencing,  this  matrix  representation 
would  be  identical  if  either  three-point  backward  or  trapezoidal  time 
differencing  were  used.  The  conditions  required  to  solve  the  blocked 
grid  problem  and  obtain  exactly  the  same  results  as  for  the  unblocked 
problem  will  be  discussed  first.  Then  some  simplifications  that  can  be 
used  and  the  error  resulting  from  these  simpl if ications  will  be 
presented.  Expanding  Equation  (U-la)  gives 

(I  ♦  -  AtRj  (4-2) 

(I  ♦  ATLn>n)Xn  -  AtRjj  -  AtLjj^Xj  (U-3) 

Equation  (4-2)  illustrates  that  the  intermediate  solution  obtained  from 
the  forward  pass  through  block  I  does  not  depend  upon  the  solution,  Xjj, 
in  block  II.  There  is  a  dependency  on  dependent  variable  values  in 
block  II,  however,  that  is  not  clearly  displayed  by  the  notation  in 
Equation  (4-2).  The  elements  of  R^  that  correspond  to  the  first  and 
second  cells  inside  the  block  boundary  are  affected  by  dependent  vari¬ 
able  values  from  time  tn  in  the  first  two  cells  across  the  block  bound¬ 
ary  in  block  II.  This  is  due  to  the  fact  that  evaluation  of  flux  at  a 
face  generally  will  require  two  points  or  either  side  of  the  face. 
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Therefore,  one  requirement  for  obtaining  exactly  the  same  answer  from 
the  blocked  calculated  as  from  the  unblocked  is  to  have  dependent  vari¬ 
able  values  available  from  the  first  two  cells  of  neighboring  blocks. 

Equation  (4-3)  shows  that  if  the  result  of  the  forward  pass  in 
block  II  is  to  be  identical  for  the  blocked  and  unblocked  cases,  then 
the  solution  from  the  forward  pass  in  block  I,  Xj ,  must  be  available. 
Though  not  shown  by  Equation  (4-3),  the  actual  elements  of  Xj  that  are 
required  are  those  in  the  first  cell  off  the  block  boundary.  Now  ex¬ 
panding  Equation  (4-lb)  for  the  backward  pass  gives 


(I  +  AtUI(I)AQj  -  X:  -  ATUIinAQIX 


d  -  ATUIIfII)AQI]:  -  xn 


(4-5) 


The  situation  here  is  similar  to  that  for  the  forward  pass,  assum¬ 
ing  that  the  correct  Xj  and  XXI  are  available.  The  backward  pass  can  be 
completed  first  in  block  II  to  obtain  AQXI  without  results  from  block  I. 
The  backward  pass  can  then  be  performed  in  block  I  to  obtain  AQj  since 
AQn  is  available. 

The  sequence  of  operations  required  to  reproduce  unblocked  results 
is  summarized  in  Figure  (4-3),  and  Figure  (4-4)  shows  how  the  discussion 
above  could  be  extended  to  more  complicated  arrangements  of  blocks.  The 
grids  depicted  by  Figures  (4-3)  and  (4-4)  can  be  considered  to  be  in  the 
computational  domain.  The  transformation  that  maps  the  blocks  back  to 
physical  space  could  map  these  rectangular  arrays  of  points  into  a  '  C 
grid  or  an  'O'  grid  or  any  other  conceivable  shape  obtainable  by  warping 
or  stretching  the  original  rectangle.  Although  it  is  possible  to  write 
the  coordinate  transformation  for  each  block  separately,  notning  would 
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be  gained  by  doing  so.  Any  blocked  grid  that  maps  into  computational 
space  in  such  a  way  as  to  allow  assembly  of  the  blocks  with  all  block- 
to-block  interfaces  correctly  aligned  can  be  represented  by  a  single 
coordinate  transformation,  except  in  the  case  where  block  faces  are  dou¬ 
ble-valued  to  allow  a  solid  surface  to  be  enclosed  between  blocks.  An 
example  where  blocking  the  grid  would  not  be  necessary  except  to  allow 
impermeable  surface  boundary  conditions  on  a  missile  fin  is  given  in 
Figure  (4-5).  For  the  class  of  problems  for  which  a  single  coordinate 
transformation  could  be  used  to  generate  the  entire  grid,  blocking  is 
still  a  useful  concept.  Computer  memory  sizes  are  not  projected  to  be¬ 
come  infinite  for  a  few  years  yet,  and  until  this  occurs,  blocked  grids 
allow  very  large  problems  to  be  solved  with  one  piece  at  a  time  in  mem¬ 
ory.  Also,  a3  discussed  in  Section  4.3,  most  charging  algorithms  heav¬ 
ily  weight  memory  usage  so  that  it  is  economically  attractive  to  reduce 
the  amount  of  memory  used  even  if  total  computational  time  increases  as 
a  result. 

The  original  motivation  for  this  study  of  blocked  grids  was  the 
capability  of  gridding  complex  geometries.  The  advantage  of  blocked 
grids  for  complex  geometries  is  the  ability  to  embed  a  block  of  one  to¬ 
pology  inside  another  block  with  a  totally  different  topology.  (Here 
the  word  topology  is  loosely  U3ed  to  refer  to  the  type  of  grid  system  as 
typified  by  the  location  of  any  cuts  necessary  to  unwrap  the  physical 
grid  into  a  rectangular  block  in  computational  space  with  no  interior 
boundaries.)  Figure  (4-1c)  is  an  example  of  this  type  of  blocned  grid 
system.  The  blocks  that  make  up  this  grid  and  their  correspondences  are 
shown  in  Figure  (4-6).  For  this  grid  and  for  most  blocked  grids  about 


truly  complex  geometries,  there  is  no  equivalent  one-block  or  unblocked 
grid.  For  many  block  arrangements  no  sequence  of  operations  exists  that 
will  provide  values  at  boundaries  when  necessary.  The  question  of 
whether  the  multiblock  solution  can  reproduce  the  unblocked  solution  is 
moot  for  these  grids,  and  the  correct  question  becomes  whether  it  is 
possible  to  transfer  information  between  blocks  in  such  a  way  as  to  pro¬ 
duce  a  stable,  conservat i ve ,  and  consistent  scheme  with  no  degradation 
in  accuracy. 

Equations  (4-ia,b)  apply  to  the  oase  where  there  exists  a  global 
unblocked  grid  equivalent  to  the  blocked  grid,  and  hence  are  not 
strictly  applicable  to  cases  where  such  an  unblocked  grid  does  not  ex¬ 
ist.  Because  blocked  grids  have  complete  continuity  across  block  bounda¬ 
ries  it  will  always  be  possible,  however,  to  trace  any  grid  line  until 
it  ends  at  a  far-field  boundary,  an  impermeable  surface,  a  reflection 
plane,  or  closes  on  itself.  Since  multiple  lines  can  originate  from  a 

single  line  as  it  passes  through  a  block  boundary this  tracing  is  not 
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unique  for  grid  lines.  As  noted  by  Norton,  Thompkins,  and  Haimes, 
however,  if  cell  centers  are  traced  then  ambiguities  are  resolved  and 
the  tracing  is  unique.  For  schemes  with  all  unknowns  along  a  computa¬ 
tional  line,  such  as  the  six-factor  and  three-factor  methods,  it  is  pos¬ 
sible  to  write  a  set  of  equations  for  each  of  the  three  lines  passing 
through  each  cell  and  solve  these  as  if  no  block  boundaries  were  pres¬ 
ent  .  This  fact  has  little  practical  value  in  terms  of  obtaining  a  solu¬ 
tion  for  the  two-factor  scheme,  since  solving  these  lines  in  the  correct 
order  to  have  necessary  information  available  for  the  solution  of  the 
next  line  ;3  net  always  possible.  It  does  mean  that  when  considering 


any  particular  block-to-block  interface  the  grid  can  be  locally  re- 
blocked  so  that  both  blocks  become  part  of  the  same  new  block.  This  is 
the  rationale  for  using  Equations  (4-ia,b)  as  a  standard  of  correct 
block-to-block  information  transfer,  even  though  an  equivalent  globally 
unblocked  grid  may  not  exist. 

As  a  very  simple  example  of  this  consider  the  one-block  "  0"  grid 
schematic  in  Figure  (4-7).  To  map  the  physical  grid  in  Figure  (4-7a) 
into  the  rectangular  block  in  the  computational  domain  shown  in  Figure 
(4-7b),  a  cut  must  be  made  along  line  ab  which  is  identical  to  cd  in  the 
physical  domain.  The  boundary  conditions  that  should  be  used  on  faces 
ab  and  cd  fall  into  the  general  class  of  block-to-block  boundary  condi¬ 
tions  under  discussion.  The  only  peculiarity  of  this  case  is  that  the 
block  on  both  sides  of  the  boundary  is  the  same.  This  simple  one  block 
example  illustrates  the  difficulty  in  supplying  correct  values  at  block 
boundaries.  In  order  to  sweep  from  face  ab  to  cd,  the  solution  at  face 
cd  must  be  used  as  a  boundary  condition  at  face  ab.  Following  the  ra¬ 
tionale  described  in  the  previous  paragraph,  the  grid  can  be  reblocked 
by  placing  the  cut  along  line  ef  instead  of  ab  making  ab  an  interior 
boundary  to  the  new  block.  If  no  error  is  incurred  on  the  outer  bound¬ 
ary  of  the  new  block  then  Figure  (4-7c)  is  identical  to  Figure  (4-2), 
and  therefore  Equation  (4-1)  is  applicable.  Of  course,  the  fallacy  in 
this  argument  is  that  boundaries  ef  and  gh  have  the  same  problem  that  ab 
and  cd  had  originally,  and  therefore  there  will  be  some  error  at  these 
boundaries.  This  is  why  the  following  arguments  concerning  error  at 
block-to-block  boundaries  can  only  be  considered  valid  in  a  small  region 
local  to  the  boundary. 
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The  value  of  Xj  in  Equation  (4-3)  and  of  AQjj  in  Equation  (4-4)  are 
the  quantities  that  will  not  necessarily  be  available  in  cases  with  no 
unblocked  equivalent.  The  simplest  approximation  for  these  quantities 
is  to  set  their  coefficients,  Uj  jj  and  Ljj  j  to  zero.  To  evaluate  the 
effect  of  this,  first  determine  the  correct  solution  for  AQj  and  AQTj. 
From  Equations  (4-2)  and  (4-3) 

X:  »  At ( I  +  AtLj  j)-1  Rj  (4-6) 

Xn  -  Ax(I  +  [Rjj  -  AtLn  j(I  ♦  AtLt^)"1  Rjj  (4-7) 

Substituting  (4-7)  in  (4-5)  and  solving  for  AQjj  gives 
AQ:I  =  (I+AtUii<ii)”1(I+AtLIIiII)~1[AtRii-At2LIIiI(I+AtL1>i)'1Ri]  (4-S) 
Using  Equations  (4-6)  and  (4-8)  in  Equation  (4-4)  gives 
AQt  -  (I  ♦  AtU^j.)"1  {(I-AtLIjI)“1AtRi  (4-9) 

-  At2Uj  t  Tj(  I+AtUji  ,  jj)-1  (I+ATLIIf  jj)-1  [RjI-AtLtji  j(  I+AtLi>  } 

Setting  Uj  jj  and  Ljj  j  to  zero  in  Equations  (4-8)  and  (4-9)  yields 

AQrr  -  At  ( I  +  Atuiitn)_1  d  +  AtLII,II)_1rII  (4-10) 
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for  the  error  in  AQjj.  To  determine  the  order  of  magnitude  of  E^Qr  and 

Eao  ,  the  magnitude  of  the  individual  factors  in  Equations  (4-12)  and 
awII 

(4-13)  must  be  known.  Each  element  of  the  matrices  Uj  j,  Uj  jj,  (Jjj  jj, 
Lj  j,  etc.  represents  one  term  of  a  first  order  spatial  difference  in 
one  of  the  three  coordinate  directions,  and  therefore  each  element  is 
either  OO/AO,  0(1/An),  or  0(1/Ac).  The  matrices  (1+AtU-j-  j)_1, 
(I+AtLj  j)  ' ,  etc.  are  0(1)  and  Rj  and  Rjj  are  0(1).  Using  this  in 
Equations  (4-12)  and  (4-13)  results  in 
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(4-1 4a) 


(4-1 4b ) 


Since  AQj  and  AQjj  are  themselves  0(At),  the  truncation  er  -or  due  to  not 
passing  the  solution  vector  across  block  boundaries  is  0(At/AE,  At/An, 
At/AO  . 

For  better  understanding  of  the  error  due  to  not  passing  the  im¬ 
plicit  solution  across  block  boundaries,  consider  the  one-dimensional 
backward  Euler  scheme,  Equation  (3-25),  rewritten  here  for  convenience. 

[I  ♦  AtA^jAQ0  -  -At«eF  ♦  0(AtAx2,  At2Ax,  At2)  (4-15) 

Taking  61  33  a  backward  spatial  difference  and  dividing  through  by  At 
gives 


AQJ  +  "J*"J 
I T  - 


A^AQ1?  -  An 


AQ, 


J-1  J-1 


Ax 


-5eFn  <•  0(Ax2,  AtAx,  At) 


(4-16) 


setting  A 


n 


J-1 


0  yields 
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((^) ’’  ♦  O(At) )  *  0(££)  -  ~((^)n+1  ♦  0(At,  Ax2)) 


2 

giving  a  truncation  error  of  0(At,  At/Axr  Ax  ) .  The  three  point  back¬ 
ward  scheme,  Equation  (3-24)  is 


3AQ"  -  AQf1  A^AQ1?  -  A?  AQ1?  . 

- J.  _ J —  ♦  J — 1 -J-J -  -5eFn  *  0(Ax2,  AtAx,  At2) 


(4-18) 


Setting  Aj_,  -  0  yields 


3aq"  -  AQj-1 


a"aq" 
+  Tx- 


-5eFn 


(4-19; 


((^)n+1  ♦  0(At2) )  +  0(£J)  -  -((H)  ♦  0(  At ,  Ax2)) 


giving  the  same  truncation  error  on  block  boundaries  as  the  backward 
Euler . 

It  is  obvious  from  the  analysis  above  that  not  transferring  the 
solution  vector  across  block  boundaries  results  in  an  inconsistent 
scheme  at  the  boundary  unless  At/Ax+0.  There  is  also  a  conservation 
error  at  the  boundary  of  0(At/Ax)  due  to  the  "implicit  flux,"  AAQ,  being 
defined  differently  in  the  adjacent  blocks.  The  conservation  error 
could  be  eliminated  by  setting  AAQ  -  0  in  both  blocks  at  the  boundary 
face,  but  this  adds  another  cell  with  an  inconsistent  differencing. 

Another  approach  for  passing  the  solution  vector  across  block 
boundaries  i3  to  approximate  Xj  in  Equation  (4-3)  and  AQjj  in  Equation 
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( 4-4 )  using  whatever  information  is  currently  available  from  other 
blocks.  The  values  used  in  this  study  are  the  results  from  both  passes 
in  adjacent  blocks,  at  whatever  time  level  was  most  recently  calculated 
in  those  blocks.  This  means  that  AQj  or  AQj-1  may  be  used  to  approxi¬ 
mate  Xj  in  Equation  (4-3),  and  A®Il  or  A®II  ^  may  be  used  to  approximate 
4air  in  Equation  (4-4),  depending  on  the  order  in  which  the  blocks  are 
solved . 

The  effect  of  using  a  time  lagged  value  of  AQ  is  easily  illustrated 
for  the  one-dimensional  backward  Euler  method.  Using  AQn-AQn_1 +0(At2) 
in  Equation  (4-16)  yields 


AQn  A?AQn  -  A"  AQ? 


j  +  j 
IF-  — 


J-l  J-1 


-<5eFn  ♦  0(Ax2,  AtAx,  At,  ¥~) 


(4-20) 


The  time  lagged  AQ  resulted  in  additional  truncation  error  of  0(At2/Ax). 

Using  the  solution  vector,  AQj,  to  approximate  the  intermediate 
solution,  Xj,  in  Equation  (4-3)  also  results  in  an  0(At2/Ax)  truncation 
error  in  the  calculation  of  Xjj  and  in  the  subsequent  calculation  of 
»«ii'  as  shown  below.  From  Equation  (4-4), 


Xj  m  AQj  *  At(Uj  jAQi  +  Uj  jjAQjj)  (4—21) 

This  matrix  equation  is  represented  in  operator  notation  by 

X  -  AQ  ♦  Ax«lAAQ  (4-22) 


which  yields 


Xj  -  AQr  ♦  0(At2) 


(4-23) 


Using  Equation  (4-23)  in  Equation  (4-3)  and  solving  for  Xjj  results  in 


XII  -  At  ( I  ♦  AtLjj^jj)  1  [*»!!  “  LII,iAQl]  +  0(T£--  JTT’ 


(4-24 ) 


8 


The  truncation  error  of  0(At2/AE,  Ai2/An,  At2/Ac)  in  Xjj  will  obviously 
carry  over  as  the  same  order  of  error  when  Equation  (4-24)  is  used  in 
Equation  (4-5)  to  calculate  AQjj.  The  order  of  the  error  is  not 
affected  if  Xn  is  approximated  by  AQn“1  since  AQn  -  AQn_1  ♦  0(At2),  and 
therefore  Equation  (4-23)  still  holds. 

To  summarize  the  results  of  the  previous  error  analysis,  setting 
the  solution  vector  values  passed  to  other  blocks  equal  to  zero  results 
in  0(At/Ax)  truncation  error,  while  approximating  the  required  solution 
vector  values  by  whatever  values  are  currently  available  from  adjoining 
blocks  results  in  0(At  /Ax)  truncation  error.  In  both  the  AQ  -  0  case 
and  the  AQ  approximated  case,  there  is  a  conservation  error  due  to 
evaluating  the  flux  at  the  block  boundary  differently  in  adjacent  cells. 

Throughout  the  above  analysis  the  assumption  was  made  that  all  val¬ 
ues  entering  the  explicit  flux  balance,  6eFn,  were  correctly  synchro¬ 
nized  to  be  time  level  n  data.  Although  it  is  not  difficult  to  code  the 
algorithm  so  that  this  is  true,  for  many  of  the  calculations  in  this 
study  the  flux  vectors  were  formed  using  the  latest  available  informa¬ 
tion  from  other  blocks.  This  technique  is  referred  to  as  using  unsyn- 
chronlzed  dependent  variables.  Using  some  time  level  n+1  information  to 
form  the  difference  5eFn,  that  should  be  totally  from  time  level  n  data 
introduces  an  0(At/Ax)  error  and  a  conservation  error  at  the  block 
boundary  face. 

For  steady  state  calculations  all  of  the  errors  discussed  above 
vanish  as  steady  state  is  approached  because  they  are  all  proportional 

to  AQ. 
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4.3  Coding  for  Blocked  Grids 

The  design  principles  used  for  writing  the  unsteady,  blocked  grid, 
implicit  Euler  code  were,  in  order  of  priority: 

1.  To  allow  any  possible  arrangements  of  blocks  with  no  change  to 
the  code. 

2.  Allow  specification  of  any  type  of  boundary  condition,  i.e.  im¬ 
permeable  surface,  far-field,  block-to-block,  etc.  on  any  patch 
of  a  block  surface. 

3.  Maximize  vectorization. 

4.  Minimize  memory  requirements. 

5.  Attain  most  efficient  use  of  high  speed  secondary  memory  (SSD  on 
CRAY  X-MP). 

One  important  advantage  of  blocked  grids  is  the  ability  to  solve 
large  oroblems  by  working  on  only  one  block  of  the  problem  at  a  time 
while  all  the  rest  of  the  data  is  stored  on  disk  or  a  similar  device.  To 
take  full  advantage  of  this,  the  code  should  have  adjustable  dimension 
arrays  so  that  arrays  associated  with  blocks  of  any  dimension  can  be 
packed  into  memory  with  no  wasted  space.  One  way  of  accomplishing  this 
in  the  FORTRAN  computer  language  is  the  ribbon  vector  method.  To  use 
the  ribbon  vector  method,  the  program  calculates  the  equivalent  one  di¬ 
mensional  array  length  for  all  adjustable  dimension  arrays.  Starting 
locations  are  then  assigned  for  each  of  these  arrays  in  the  blank  common 
block,  making  sure  to  allow  enough  space  between  the  start  of  one  array 
and  the  next  to  store  the  entire  array.  Most  computer  systems  provide  a 


user  callable  subroutine  that  will  allow  the  user  to  set  the  amount  of 
memory  available  to  the  program.  This  process  is  referred  to  here  as 
dynamic  memory  management.  On  systems  with  dynamic  memory  management 
blank  common  should  contain  one  array  singly  dimensioned  to  a  length  of 
one.  After  the  user  program  has  determined  the  total  number  of  words 
required  in  the  ribbon  vector,  the  program  will  call  the  dynamic  memory 
management  routine  to  allocate  that  many  extra  words.  Because  blank 
common  is  stored  as  the  last  item  of  the  users  program  in  memory,  the 
extra  allocated  space  goes  to  blank  common. 

Specific  arrays,  singly  or  multiply  dimensioned,  in  the  ribbon  vec¬ 
tor  can  be  easily  referred  to  and  used  inside  subroutines  by  passing  the 
array  as  a  parameter  in  the  calling  statement  and  using  standard  FORTRAN 
adjustable  dimension  arrays  in  the  subroutines.  Figure  (4-8)  gives  an 
example  of  ribbon  vector  usage. 

As  discussed  in  Section  4.2,  certain  Information  must  be  trans¬ 
ferred  between  blocks  to  obtain  reasonable  results.  This  information  is 
invariably  in  the  first  two  cells  of  neighboring  blocks  and  may  consist 
of  the  dependent  variable  values  only  or  the  dependent  variables  and  the 
most  recently  available  solution  vector  value.  Since  only  one  block  is 
contained  in  memory  at  a  time,  special  provisions  must  be  made  for 
transferring  needed  values  from  other  blocks.  The  approach  used  here  is 
to  store  a  set  of  resident  boundary  condition  arrays  permanently  in  mem¬ 
ory  for  each  block  surface  patch  that  joins  another  block.  These  resi¬ 
dent  arrays  are  updated  to  the  next  time  level  from  information  in  each 
block  while  that  block  is  in  memory.  By  doing  this  it  is  never  neces¬ 
sary  to  read  or  write  a  block  more  than  once  per  iteration  through  the 


entire  grid.  Each  block  is  dimensioned  to  have  enough  phantom  points  to 
store  all  boundary  condition  information.  Block  boundary  condition  in¬ 
formation  is  copied  into  the  block's  phantom  point  locations  from  the 
resident  boundary  condition  arrays  for  subsequent  calculations.  This 
greatly  simplifies  vectorization  of  calculations  as  opposed  to  using 
logic  to  pull  points  from  the  permanent  boundary  arrays  as  needed. 

Doubly  storing  the  dependent  variable  values  in  this  way  also  al¬ 
lows  synchronization  of  terms  used  in  the  explicit  flux  balance.  The 
method  used  to  insure  that  time  level  n  data  is  always  used  at  block 
boundaries  is  given  below.  As  a  starting  condition  let  all  blocks  and 
the  resident  BC  (Boundary  Condition)  arrays  contain  level  n  data.  The 
algorithm  is  as  follows: 

1 .  Bring  current  block  into  memory  and  update  block  BC  phantom 
points  from  all  resident  BC  arrays. 

2.  Store  information  from  current  block  into  all  resident  BC  arrays 
associated  with  higher  numbered  blocks. 

3.  Solve  for  Qn+1  in  current  block. 

M.  Store  information  from  current  block  into  all  resident  BC  arrays 
associated  with  blocks  numbered  less  than  the  current  block  and 
resident  BC  arrays  for  the  current  block. 

The  high  speed  secondary  memory  on  the  CRAY  X-MP,  referred  to  as 
SSD,  is  used  to  store  all  blocks  not  currently  in  memory.  The  SSD  is 
accessed  using  standard  FORTRAN  input  and  output  methods.  To  make  the 
most  efficient  use  of  SSD,  all  SSD  data  files  were  unblocked,  BUFFER  IN 


and  BUFFER  OUT  statements  were  used  to  transfer  all  data  associated  with 
a  block  in  a  single  unformatted  buffer  command. 

4.4  Two-Dimensional  Calculations  on  Blocked  Grids 

This  section  presents  results  from  steady  and  unsteady  calculations 
on  blocked  grids.  Convergence  and  stability  characteristics  are  exam¬ 
ined  for  the  steady  case,  and  the  converged  steady  solution  is  shown  to 
not  depend  on  the  blocking  of  the  grid.  The  error  in  unsteady  results 
due  to  block  boundary  positioning  is  studied  for  the  same  airfoil  and 
flow  conditions  presented  in  Section  3-5  and  results  with  AQn-0  at  block 
boundaries  and  also  with  AQn  approximated  by  whatever  information  is 
available  from  adjoining  blocks  are  presented. 

The  steady  calculations  were  for  a  NACA0012  airfoil  at  Mach  0.80 
and  angle  of  attack  of  1.25  degrees.  The  grid  used  first,  and  the  grid 
that  would  later  be  split  up  to  form  block  grids  for  comparison,  was  the 
one-block  "C"  grid  described  in  Section  3-5  for  use  in  unsteady  calcula¬ 
tions.  Figure  (4-9)  shows  the  pressure  distribution  obtained  after  1000 
local  time  steps  at  a  CFL  of  15  using  backward  Euler  time  differencing. 
The  strong  upper  surface  shock  was  sharply  captured,  but  the  lower  sur¬ 
face  3hock  i3  not  properly  resolved.  Clustering  grid  points  in  the 
neighborhood  of  the  expected  lower  surface  shock  would  sharpen  it,  but 
this  was  felt  to  be  unnecessary  for  this  study  of  blocked  versus  un¬ 
blocked  results.  The  blocked  grid  tests  were  chosen  to  provide  worst 
cases  by  placing  block  boundaries  in  regions  where  the  flowfield  is  rap¬ 
idly  changing.  Figure  (4-10)  illustrates  the  division  of  the  original 
grid  into  four  blocks  by  cutting  the  grid  at  the  leading  and  trailing 


edges.  This  is  a  critical  positioning  of  block  boundaries  due  to  the 
high  gradients  in  the  leading  and  trailing  edge  stagnation  regions,  par¬ 
ticularly  during  startup  from  the  initial  conditions  of  free-stream  val¬ 
ues  everywhere.  The  first  block  boundary  treatment  used  was  to  set  AQn 
»  0  at  the  boundary  and  to  use  the  latest  available  (unsynchronized) 
dependent  variable  values  from  adjoining  blocks.  Although  the  starting 
transients  for  this  case  produced  highly  unrealistic  pressure  values 
particularly  at  the  leading  edge,  the  calculation  remained  stable  using 
local  time  stepping  at  a  CFL  of  15,  and  steady  results  identical  to  the 
one-block  solution  were  obtained.  Figure  (4-ii)  shows  a  three-block 
grid  chosen  to  investigate  possible  interaction  between  shock  waves  and 
the  block  boundaries.  Steady  results  for  this  grid  were  also  identical 
to  the  one-block  case  using  the  same  block  boundary  treatment.  These 
tests  show  that  stability  is  insensitive  to  block  boundary  positioning 
and  that  no  error  is  introduced  in  the  steady  solution  due  to  blocking. 

The  convergence  history  for  the  one-block,  three-block,  and  four- 
block  steady  calculations  as  measured  by  the  maximum  change  in  density 
is  plotted  in  Figure  (4-12).  Another  indicator  of  convergence  is  the 
iteration  number  at  which  the  number  of  supersonic  points  stops  chang¬ 
ing,  or  is  frozen.  The  number  of  supersonic  points  for  the  one-block 
case  is  frozen  at  474  Iterations;  for  the  four-block  case  at  604  itera¬ 
tions;  and  for  the  three-block  case  at  599  iterations.  For  the  cases 
considered  approximately  25J  more  iterations  are  required  to  obtain  the 
same  degree  of  convergence  with  the  blocked  grids  as  with  an  unblocked 
grid.  This  disadvantage  is  more  than  offset  in  terras  of  cost  to  obtain 
a  solution,  however,  by  the  greatly  reduced  memory  requirements  for 
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blocked  grid  calculations.  Table  1  shows  a  comparison  of  various  job 
statistics  from  the  three  diffferent  cases  run  on  the  NASA-Ames  CRAY 
X-MP.  Note  that  even  though  CPU  times  for  the  blocked  grid  cases  were 
higher  than  the  unblocked  cases,  the  cost  in  terms  of  accounting  units 
per  step  for  the  blocked  grid  cases  were  one-third  that  of  the 
unblocked. 

Instead  of  setting  AQ  -  0  and  using  the  unsynchronized  dependent 
variable  values  at  block  boundaries,  the  solution  vector,  AQ,  could  be 
approximated  and  synchronized  dependent  variable  values  could  be  used. 
In  Section  4.2  it  was  shown  that  this  should  produce  less  error  at  block 
boundaries  for  unsteady  problems.  The  effect  on  convergence  to  steady 
state  of  this  more  accurate  boundary  treatment  is  illustrated  in  Figure 
(4-13)  for  the  three-block  case.  The  combination  of  unsynchronized  de¬ 
pendent  variables  with  AQ  approximated  is  not  shown  because  it  was  un¬ 
stable  at  a  CFL  of  15.  The  more  accurate  boundary  treatment  obviously 
does  not  improve  convergence  to  steady  3tate.  Its  effect  on  time-accu¬ 
racy  will  be  discussed  below. 

The  same  foui — block  and  three-block  NACA0012  grids  used  for  the 
steady  calculations  discussed  above  were  used  to  calculate  unsteady  flow 
at  Mach  0.755  and  a  reduced  frequency  of  0.1628.  The  mean  angle  of  at¬ 
tack  was  0.016  degrees  and  the  amplitude  of  the  unsteady  pitching  motion 
was  2.51  degrees.  These  are  the  same  conditions  used  in  Section  3.5. 
The  time  step  was  taken  to  be  At  -  0.05,  and  backward  Euler  time 
differencing  was  used  for  all  unsteady  block  calculations. 

Pressure  distribution  for  the  one-block,  three-block  and  four-block 
grids  are  shown  in  Figures  (4-1 4a)  through  (4-i4h).  These  results  were 
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obtained  by  using  unsynchronized  dependent  variables  from  adjoining 
blocks  and  setting  AQ  -  0  at  the  block  boundary.  The  four  block  results 
are  nearly  indistinguishable  from  the  one  block,  but  the  three-block 
results  show  a  considerable  error  in  shock  location  everywhere  except  in 
Figures  (4-14d)  and  (4-14h).  In  Figures  (4-l4d,h)  the  shock  wave  has 
reached  its  position  of  maximum  downstream  travel  on  the  upper  and  lower 
surface,  respectively,  and  is  nearly  stationary.  The  shock  wave  posi¬ 
tion  given  by  the  three-block  calculation  lags  the  correct  location  in 
all  cases  that  have  an  appreciable  shock  velocity.  The  excellent  agree¬ 
ment  of  the  four  block  case  is  due  the  small  variation  with  time  depend¬ 
ent  variables  at  the  block  boundaries.  Even  though  large  spatial 
gradients  are  present  at  the  leading  and  trailing  edge  of  the  airfoil, 
the  pitching  motion  does  not  cause  rapid  changes  in  a  particular  cell  in 
this  region. 

The  analysis  of  Section  (4.3)  showed  that  by  approximating  AQ  at 
block  boundaries  and  synchronizing  the  dependent  variables  used  from 
adjoining  blocks,  the  error  at  the  boundary  could  be  reduced  from 
0(At/Ax)  to  0(At2/Ax).  Figures  ( U — 1 5a )  through  (4-l5h)  show  a  compari¬ 
son  between  pressure  distributions  obtained  using  approximated  AQ'3  with 
synchronized  dependent  variables,  AQ-0  with  synchronized  dependent  vari¬ 
ables,  and  AQ-0  at  the  boundary  and  unsynchronized  dependent  variables. 
The  combination  of  approximated  AQ  with  unsynchronized  dependent  vari¬ 
ables  was  not  obtained  due  to  stability  problems.  The  AQ  -0  with  unsyn¬ 
chronized  dependent  variables  results  are  the  same  as  those  in  Figure 
(4-14).  Using  synchronized  dependent  variables,  but  maintaining  the 
rough  AQ-0  approximation  did  not  significantly  improve  the  comparison 


with  one-block  results.  For  example,  at  the  instant  in  time  shown  in 
Figure  (4-1 5a),  the  AQ-0  with  synchronized  dependent  variables  compari¬ 
son  is  not  as  good  as  the  AQ-0  with  unsynchronized  dependent  variables. 
Examination  of  Figure  (4-15)  reveals  the  dramatic  improvement  yielded  by 
the  approximated  AQ  and  synchronized  dependent  variable  block  boundary 
conditions.  The  largest  difference  between  the  one  block  results  and 
results  from  using  approximated  AQ  and  synchronized  dependent  variables 
occurs  when  the  lower  surface  shock  moves  upstream  from  block  one 
through  the  block  boundary  roughly  parallel  to  the  3hock  front  and  into 
blocks  two  and  three  (See  Figure  ( 4— 1 1 )  for  the  grid  geometryj.  The 
lower  surface  shock  has  Just  left  block  one  in  Figure  (4-15a).  The  up¬ 
per  surface  shock  in  Figure  (4-1 5e)  has  roughly  the  same  relative  posi¬ 
tion  and  motion  as  the  lower  surface  shock  in  Figure  (4-l5a),  but  for 
this  case  the  approximated  AQ  results  show  no  error  in  shock  position. 

The  results  presented  in  this  section  from  steady  and  unsteady  cal¬ 
culations  on  blocked  grids  have  demonstrated  the  effects  on  convergence 
to  steady  state  and  on  time-accuracy  of  various  block  boundary  condi¬ 
tions.  All  of  the  block  boundary  conditions  represented  led  to  the  same 
steady  solution,  but  the  convergence  rate  was  slower  on  blocked  grids 
than  unblocked  grids.  Block  boundary  conditions  that  are  formally  more 
time-accurate  slowed  convergence  more  than  the  less  accurate  conditions 
for  the  cases  presented  here.  Cost  comparisons  for  the  steady  calcula¬ 
tions  showed  that  the  reduced  memory  requirements  of  blocked  grids  make 


their  use  attractive  even  if  more  computer  time  must  be  spent  to  obtain 
the  same  level  of  convergence.  Comparison  of  unsteady  calculations  on 
blocked  and  unblocked  grids  verified  earlier  analysis  that  showed  that 
approximating  AQ  on  the  block  boundaries  and  using  synchronized  values 
of  the  dependent  variables  introduces  less  error  than  setting  AQ-0  and 
using  unsynchronized  dependent  variables. 


CHAPTER  V 


THREE-DIMENSIONAL  CALCULATIONS 


This  chapter  presents  a  comparison  of  steady  and  unsteady  calcula¬ 
tions  with  experimental  results  for  a  rectangular  supercritical  wing  of 
aspect  ratio  two.2^  The  grid  for  the  wing  was  generated  as  a  single 
continuous  block  of  100x27x25  points  and  then  split  into  four  blocks  as 
shown  in  Figure  (5-1).  The  overall  grid  is  a  "C"  mesh  in  the  streamwise 
direction  and  half  an  "0"  mesh  in  the  spanwise  direction.  Points  are 
clustered  near  the  leading  edge  and  near  the  wing  tip.  Blocks  1  and  3 
include  most  of  the  lower  and  upper  surfaces  of  the  wing,  respectively, 
while  block  2  wraps  around  the  wing  tip.  Each  of  these  three  blocks  are 
81x27x9.  Block  4  includes  all  points  downstream  of  the  wing  and  is 
20x27x25.  The  wind  tunnel  test2^  was  conducted  in  freon  so  that  Y  - 
1.131  for  the  calculations. 

Figure  (5-2)  shows  steady  state  pressure  results  at  four  span  sta¬ 
tions  for  which  experimental  data  were  available  for  a  Mach  number  of 
0.701  and  four  degrees  angle  of  attack.  The  semispan  location  given  on 
the  figures  in  this  chapter  is  for  the  center  of  the  cell  used  to  com¬ 
pare  with  experiment.  Even  though  the  span  location  of  the  experimental 
data  does  not  exactly  match  the  cell  center  location,  the  experimental 
data  point  is  well  within  the  boundaries  of  the  computational  cell  with 
which  it  is  compared  in  all  cases.  The  calculated  leading  edge  expan¬ 
sion  overshoots  the  experimental  near  the  wing  root,  but  slightly  under¬ 
expands  at  the  tip.  Agreement  at  the  tip  is  better  than  at  the  root, 
but  the  inviscid  calculation  generally  predicts  lower  upper  surface 
pressures  than  were  measured.  Steady  c<.  ‘ culations  were  done  for  this 
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Mach  number  at  two  degrees  angle  of  attack  with  similar  results.  At 
higher  Mach  numbers  the  comparison  with  experimental  results  worsens. 
For  Mach  0.83  and  two  degrees  angle  of  attack  the  Euler  code  predicts  a 
shock  wave  at  the  trailing  edge  of  the  wing,  but  the  experiment  shows 
the  shock  much  farther  forward.  Figure  (5-3)  is  a  plot  of  experimental 
and  calculated  results  for  this  case  at  the  60$  semispan  location.  The 
fairly  strong  three-dimensional  effects  present  on  this  low  aspect  ratio 
wing  are  modelled  fairly  well,  but  viscous  effects  need  to  be  accounted 
for  to  significantly  improve  agreement  with  experiment. 

Unsteady  calculations  for  oscillatory  rigid  body  pitching  of  the 
wing  were  performed  to  compare  with  experimental  results.  The  pitch 
axis  of  the  wing  was  located  at  46$  of  the  chord,  and  the  amplitude  of 
the  unsteady  angle  of  attack,  a, ,  was  one  degree.  Two  reduced  frequen¬ 
cies,  k-0.358  and  k-0.714,  were  studied  at  a  freestream  Mach  number  of 
0.7  and  a  mean  angle  of  attack,  aQ,  of  four  degrees.  The  time  step  size 
used  for  both  reduced  frequencies  was  approximately  At  -  0.07.  This  time 
step  size  resulted  in  a  maximum  CFL  of  approximately  130  in  blocks  1  and 
3  and  approximately  60  in  block  4.  In  the  wing  tip  block,  block  2,  the 
maximum  CFL  was  around  500  due  to  the  very  small  cells  near  the  line 
singularity  leaving  the  leading  and  trailing  edge.  The  time  step  sizes 
chosen  resulted  in  360  time  steps  per  cycle  of  motion  for  k-0.358  and 
180  time  steps  per  cycle  of  motion  for  k-0.714.  The  magnitude  and  phase 
of  the  unsteady  pressure  coefficients  were  obtained  by  Fourier  analysis 
of  the  last  cycle  of  motion  calculated  at  each  frequency.  Figure  (5-4) 
presents  magnitude  and  phase  results  obtained  from  the  third  cycle  of 
motion  at  k-0.358.  Two  hundred  local  time  step  iterations  at  a  CFL  of  8 


preceded  the  start  of  the  unsteady  motion.  The  flow-field  information 
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from  the  last  time  step  of  this  calculation  was  used  to  begin  the  calcu¬ 
lation  at  k-0.714,  and  the  higher  frequency  motion  was  continued  for 
four  cycles  of  motion.  Magnitude  and  phase  of  the  unsteady  pressures  at 
k— 0 - 7 1 ^  are  plotted  in  Figure  (5-5). 

Figures  (5-4)  and  (5-5)  show  the  characteristic  spike  in  the  magni¬ 
tude  of  unsteady  pressure  due  to  shock  wave  motion.  The  position  of  the 
computed  spike  is  different  from  the  experiment  by  the  same  amount  that 
the  placement  of  the  steady  shock  was  in  error,  as  seen  by  referring 
back  to  Figure  (5-2).  The  difference  in  amplitude  between  Euler  results 
and  experiment  is  difficult  to  determine  due  to  the  sparseness  of 
experimental  data  points,  but  the  Euler  code  seems  to  overpredict  the 
magnitude  near  the  root  and  underpredict  near  the  tip.  The  agreement  in 
phase  is  generally  very  good  except  for  the  effects  of  the  mislocated 
shock.  An  exception  to  this  is  the  82*  semispan  results,  for  which 
there  may  be  problems  with  the  experimental  results. 

Due  to  the  large  effect  on  unsteady  forces  and  moments  that  can  be 
expected  from  the  shock  spike,  a  relatively  small  error  in  shock  loca¬ 
tion  can  result  in  relatively  large  differences  in  unsteady  forces  and 
moments.  The  shock  positioning  effects  of  viscosity  are,  therefore, 
very  important  for  unsteady  transonic  analysis. 
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CHAPTER  VI 


CONCLUSIONS 


An  unsteady  implicit  Euler  equation  solution  algorithm  using  a  fi¬ 


nite  volume  discretization  and  flux-vector  splitting  has  been  presented. 


The  effect  on  time-accuracy  of  different  time  step  sizes,  different  ap¬ 
proximate  factorizations,  and  formal  first-order  versus  second-order 
time  accuracy  has  been  determined  by  numerical  experimentation  on  an 
oscillating  NACA0012  airfoil  in  transonic  flow.  It  was  found  that  time 
step  sizes  corresponding  to  Courant  numbers  of  100  or  more  can  produce 
time-accurate  results  if  flow  variables  are  not  rapidly  changing.  Due 
to  better  stability  properties,  the  two-factor  method  studied  here  gave 
better  results  than  the  six-factor  method.  The  second-order  time-accu¬ 
rate  three  point  backward  scheme  yielded  only  slight  improvement  for  the 


M 


test  cases  over  backward  Euler  time  differencing,  which  is  only 
first-order  in  time. 

Another  topic  covered  here  was  the  use  of  blocked  grids  to  obtain 
time-accurate  Euler  solutions.  Techniques  for  transfer  of  information 


between  blocks  were  analyzed  and  verified  by  comparing  multi-block  solu¬ 
tions  with  equivalent  one-block  solutions.  It  was  shown  that  approxi¬ 


mating  the  value  of  the  solution  vector  required  at  block  boundaries 
with  whatever  information  is  currently  available  from  adjoining  blocks 
introduces  an  0(At  /Ax)  error  in  the  transient  solution  at  the  boundary, 
and  gives  unsteady  results  that  compare  well  with  unblocked  results  even 
for  cases  with  a  shock  wave  passing  through  the  block  boundary. 
Converged  steady  results,  with  no  body  motion,  showed  no  error  due  to 
the  block  boundaries.  Convergence  to  steady  state  was  slowed  to  some 
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extent  by  using  blocked  grids,  but  the  use  of  blocked  grids  is  economi¬ 
cally  attractive  due  to  the  greatly  reduced  computer  memory  cost. 

Steady  and  unsteady  calculations  of  transonic  flow  over  a  rectangu¬ 
lar  planform  supercritical  wing  using  a  blocked  grid  were  compared  to 
experimental  results.  The  steady  comparison  is  good  for  those  cases  in 
which  viscous  effects  are  negligible,  but  for  higher  Mach  number  cases 
the  inviscid  calculation  places  the  shock  downstream  of  the  experimental 
location.  Unsteady  calculations  showed  good  agreement  in  terms  of  mag¬ 
nitude  and  phase  of  the  pressure  on  the  wing  except  where  the  shock  wave 


was  misplaced. 
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(a)  25°  or  Oscillatory  Motion,  o  Increasing  Throuoh  1.09° 

Figure  (3-3)  Unsteady  Airfoil  Pressure  Coefficients  Using  Various  Time  Step 
Sizes  and  fime  Discretizations 
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SHOCK  LOCATIONS 


Unsteady  Airfoil  Lift  and  Moment 
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DT1II  N=0 .10,  2nd  ORDER 
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PROGRAM  RIBBON 
COMMON  S(l) 


This  program  will  read  three  3-D  arrays  with  dimensions  NI,NJ,NR. 

READ  (5,*)  NI,  NJ,  NK 
I  SIZE  -  NI*NJ*NK 


t 

v.v.v. 

v  v 

VaV.' 

.••’A'.V 

a*  < 


LX,  LY,  and  LZ  are  the  starting  locations  for  X,  Y,  and  Z. 
LX-1 

LY-LX+ISIZE 

LZ-LY+ISIZE 

LAST"LZ-*-ISIZE-l 

Subroutine  GETMEM  will  increase  the  available  memory 
CALL  GETMEM  (LAST) 

CALL  READXYZ(S(LX)  ,S(LY ) ,  S  (LZ) ,NI , NJ,NK) 

END 

SUBROUTINE  RFADXYZ(X,Y,Z,NI,NJ,NK) 

DIMENSION  X(NI,NJ,NK) , Y (NI,NJ,NK) ,Z(NI ,NJ ,NK) 

READ(10)X,Y,Z 

DO  20  K-l.NK 

DO  20  J-l.NJ 

DO  20  1*1, NI 

A  -  X(I,J,K)  Y( I , J, K) 


20  Z( I, J,K)  *  A*Z(I,J,K) 

RETURN 

END 


Figure  (4-8)  Use  of  Adjustable  Arrays  in  FORTRAN 
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CONVERGENCE  HISTORIES 

Local  Time  Stepping  with  CPL=15 
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Figure  (*-15)  Unsteady  Pressure  Coefficients  Using  Various  Block  Boundary 
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a)  continued 


i)  continued 
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APPENDIX  A 


TIME-DEPENDENT  CURVILINEAR  COORDINATES 


Rewriting  the  nondimens ional  form  of  Equation  (2-1)  using  the 
transformed  coordinates  given  by  Equation  (2-3)  and  the  chain  rule 
yields 
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Then 
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(A-2) 


where 
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where  V*  must  depend  on  t  in  order  for  the  unsteady  coordinate  represen¬ 
tation  to  refer  to  the  fixed  volume  V g  in  Cartesian  space  and  J  is  given 
by  Equation  (A-3).  Using  Leibnitz'  rule  to  take  the  time  derivative  of 
the  integral  results  in 
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( A— 7 ) 


where  2  -  Ut ,  nt,  ?t) 

and  fi  is  the  unit  normal  in  (5tnfc)  space  to  the  surface  element  dS  of 
volume  V  (t).  Application  of  Gauss’  Theorem  to  the  second  integral 
yields 


V  ’  ( t ) 


3(JCJ 


3(Jnt) 

— + 


3(J;t). 

-y—  jdgdndc 


0. 


( A— 8 ) 


Since  the  original  volume,  VQ,  is  arbitrary,  Equation  ( A— 8 )  must  hold 
for  all  V '  ( t ) .  This  implies  that  the  integrand  must  be  identically 
zero,  yielding 
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The  first  term  in  Equation  (A-9)  can  be  rewritten  as 
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because  t  is  a  function  of  t  only.  Multiplying  through  by  tT  in  Equa¬ 
tion  (A-9)  and  using  Equation  ( A— 2 )  and  Equation  (A-10)  yields 
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Therefore,  the  coefficient  of  q  in  Equation  (A-4)  is  identically  aero. 
3ecause  Equation  (A-9),  or  equivalently,  Equation  ( A— 11),  express  con¬ 
servation  of  volume,  they  are  referred  to  as  statements  of  the  geometric 
conservation  law. 

Evaluation  of  the  coefficients  of  f,  g,  and  h  on  the  right  hand 

t  1  I 

side  of  Equation  (A-4)  requires  writing  the  quantities  J  £x.  J  nx,  J  Cx. 
etc.  in  terms  of  x^,  x^,  Xj.,  y^,  etc.  To  do  this,  first  use  the  chain 
rule  to  express  partial  derivatives  with  respect  to  the  curvilinear  co¬ 
ordinates  in  terms  of  partial  derivatives  with  respect  to  the  Cartesian 
coordinates  as 
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where  a,2  -  -[xT  (ynVycV  -  *  zT  (x^-x,/^  ] , 


a1 2  “  xt^y5z!;"'y!;z5)  “  yT  ( x^Zj.-x^z^ )  +  z^x^-x^) , 
a13  -  -CxT(y5zn-ynz?)  -  yT(x^zT)-xnz^)  +  zT  (x^-x^)  ] . 


Comparing  Equation  ( A— 1 3 )  with  the  chain  rule  expression  for  the  partial 
derivatives  with  respect  to  the  Cartesian  coordinates  in  terms  of  par — 
tial  derivatives  with  respect  to  the  curvilinear  coordinates, 


gives  the  values  of  the  metrics. 
Specifically,  the  metrics  are 
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<zztK 

( A-1 5n 1 

The  coefficients  of  f,  g,  and  h  in  Equation  (A-4)  can  now  be  evalu¬ 
ated.  Using  Equation  (A-15),  the  coefficient  of  f  can  be  written  as 


^(ynzC  “  +  In' 


lT^(y 

This  is  of  the  form 


<Ve  “  Hz0  *  lc  (y5zn  "  ynVl-  (A"16) 


tT  7  •  (7y  X  7z) , 


(A-17) 


where  the  vector  operator  7  is  taken  to  be 


7  =  f3  3  3  ] 
v  ~  l3T  7n*  3?J 


(A-18) 


and  the  •  and  X  represent  vector  dot  and  cross  products,  respectively. 
Using  a  well-known  vector  identity.  Expression  (A-17)  can  be  rewritten 
as 

tT[7z  •  (7  X  7y )  -  Vy  •  (7  X  7z)]  =  0. 


( A— 1 9 ) 
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The  expression  is  identically  equal  to  zero  because  the  curl  of  the  gra¬ 
dient  is  always  zero. 

Similarly,  the  coefficient  of  g  on  the  right  hand  side  of  Equation 
(A-4)  can  be  written  as 


tT  V  •  (Vz  X  7x5  =  0, 


(A-20) 


■  .A  j 


*  «»■  ■> * 

-f.vV;. 
.v  .y 


and  the  coefficient  of  h  can  be  written  as 


tT  V  •  (7x  X  Vy)  =  0. 


(A-21 ) 


«*  / 

•  v  * V 

s  *,  *,  • 
A  / 


Using  the  metric  identities  given  by  Equations  ( A— 1  1 )  ,  ( A— 19), 

(A-20),  and  (A-21)  in  Equation  (A-4)  yeilds 


3Q  3F  3G  4  3H  . 

TT  3n  3c  “ 


(A-22) 


r  v  V*\- 
■  .  .« 


where  Q  »  Jq 


F  -  JtT(Etq  ♦  Exf  *  Cyg  *  5zh) 
G  •  JtT(ntq  ♦  nxf  *  nyg  ♦  nzh) 
H  -  JtT(qtq  ♦  cxf  *  Cyg  *  Czh). 


»  ^ 
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APPENDIX  B 


JACOBIANS  OF  THE  SPLTT  FLUX  VECTORS 


The  Jacobians  of  the  split  flux  vectors  with  respect  to  the  depend¬ 
ent  variables  are  derived  in  this  appendix  in  a  form  suitable  for 

numerical  calculation.  The  split  flux  vectors  referred  to  in  Equation 
(3-7)  are  defined  for  the  one-dimensional  case  in  Equation  (3-31).  For 
the  three-dimensional  case  the  general  flux  vector  is 

Kl  -  l  *  IxJ(QL)  hK£(QL)  (B-1  a ) 

1-1 ,4,5  K 

KR  -  l  i(xJ(QR)  -  I xf (QR) 1 )K,(QR )  (B-1 b) 

1-1, 4, 5 

where  the  symbol  K  should  be  replaced  by  F,  G,  or  H  when  k  is  E,  n,  or 
; ,  respectively,  and  and  Q  are  obtained  by  dependent  variable 

extrapolation  along  the  k  coordinate  from  more  negative  coordinate  val¬ 
ues  and  more  positive  coordinate  values,  respectively.  The  eigenvalues, 
X,4,  and  xjj  are  given  by  Equation  (3-4),  and  the  split  flux  vectors, 

<1  ,  Ky ,  and  are  given  by  Equation  (3-6).  Equation  (B-1)  illustrates 

that  and  KR  are  the  3um  of  from  zero  to  three  components,  depending 

on  the  signs  of  the  eigenvalues. 

The  Jacobian  of  each  of  the  possible  components, 

9(AtKl) 

h  "  — — -  1  “  1'4'5*  (B~2) 


will  be  derived. 

Expanding  the  derivative  in  Equation  (B-2)  yields 
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3K1 
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The  elements  of  the  dependent  variable  vector  are  given  by 


Q1  ■  Jp 
Q2  »  Jpu 
»  Jpv 
Q4  »  Jpw 
Q5  *  Je 

From  Equation  (3-4),  the  eigenvalues  are 

*k  *  Vk 

xk  “  tt(3k  +  c^7k^  “  M9*  +  c^7k^  +  tTkt 
A-5  -  tT(0k  -  c  I  Vk  1  )  -  tT  ( 9k  -  c)  1 7k  I  *  tTkt 

where 

Sk  -  kxu  ♦  ky v  *  kzw  ♦  kt 

In,  I  ,,  2  ,2  ,  2,1/2 

iVkl  -  (kx  ♦  ky  +  kz) 

9k  *  kxu  *  kyv  +  kzw 
9k 

9k  "  -[Wf 


The  local  speed  of  sound,  c,  is 


c 


where 


Expressing  Sk,  0.K,  c,  and  p  in  terms  of  the  dependent  variables  yields 


kxQ2  ♦  kyQ3  +  kzQ4 


(B-1 1 : 


KxQ2  ♦  ky^3  ♦  kZQU 


(B-1 2) 


,  Qc  y/v.n  Q?  *  Qf  +  Qw  1/2 

r vcv  _  ii  5  _  '  1 '  2  3  4-| 

ty{y  -  1)3-  —5 - - \ 


(B-1 3) 


Q2  +  Q2  4-  q2 


(B-1 4 ) 


Note  that 


I  I  ^®k 

l7ki  TCJ*  J  "  ,-2«3.^5. 


(B-15) 


Using  Equations  (B-5),  (B-11),  (B— 12) ,  (B— 1 3) »  and  (3-15),  the 

derivatives  of  the  eigenvalues  with  respect  to  the  dependent  variables 
can  be  written  as 


9l,  k„  k„  k„ 


k  :  k  *x  *y  *z  1 

r  M~  37*  57*  57’  37*  °J 


(B-16) 


W  ”  tT^7k^'3U7(®k+c)'  nj^(Vc)’  ^(Vo)’  -J57(®k*c)*  (B-17) 

8 15 

nr  "  tT'Vk^'^7(®k~c)>  ^(®k'c)>  TO7(®k”o)»  3^(®k”c)»  ~  lijrJ  (B-1 8 ) 

The  derivatives  of  9^  are  given  by 


38.  9, 

k  k 

^7  "  "  o7 


(B-1 9a) 


39.  a 
<  x 


(B-1 9b ) 


(B-I9c) 


(B-I9d) 


( B— 1 9e ) 


where 


‘y  "  JW[ 


The  derivatives  of  c  are  most  easily  evaluated  by  rewriting  c  as  ct^/Q, 
Using  Equation  (B-13),  the  derivatives  of  cQ1 ,  are 


aCcQ^  YOr-i)Q5 


9(cQ1 )  T(V1)Q2 


aecQ^  Y(y-1)Q3 


aCcq^  y(y-i)Q4 


3(cV  t(y-i  )Q1  Y(Y_n 


(B-20a) 


(B-20b) 


(B-20C) 


(B-20d) 


( 3-20e ) 
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Tne  derivatives  of  c  are 


3c 

*37’ 

S1 

3(cQ, 

3c 

0-1 

3(cQ1  ) 

“ 

Q1 

3c 

n-1 

3(05,) 

3dj 

Q1 

3c 

n-1 

3(cQ.  ) 

3Qy 

Q1 

3c 

n-1 

3(cQ,  ) 

*5“ 

Q1 

(B-21a) 

(B-21 b ) 

(3-21c) 

(B-21 d) 

(B-21 e ) 


The  derivatives  of  pJ  with  respect  to  the  dependent  variables  will  be 
needed  below.  Using  Equation  (8-1 H),  these  values  are 


3  (pJ ) 

Y-1  Q2  *  Q3  *  «« 

(B-22a) 

i 

of 

3  (pJ) 

a 

l 

1 

(B-22b) 

3  ( pJ ) 

Q3 

'  '  37 

(B-22c) 

3(pJ) 

a 
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i 

(B-22d) 

3(pJ) 

“TOT 

-  Y-1  . 

(B-22e ) 

K,  ,  Kn ,  and  Kc  written  in  ter*ns  of  the  dependent  variables  are 


g2  +  g2  +  g2 
w2  3  4 


Q->  +  cQ,  k 


(3-2 3a) 


^  ’  77 


Q3  -  CQ^y 


(3-23b) 


Qh  * 


Q5  *  pJ  *  cQ1 9^ 


Qo  -  cQik 


'5  "  77 


Q3  -  CQ^y 


(B-23c) 


Qk  “  cQ, k. 


Qc  ♦  pJ  -  cQ,9, 


The  derivatives  of  the  split  flux  vectors  with  respect  to  the  dependent 
variables  are 


W  77  x 


3(cQ1  ) 


8(cQ1  ) 


x  “to;  x  “3$;  * 


_  3(cQ1  )  _  3(cQ1 )  _  3(cQ  ) 

'*y  a02  1_^  ~ky  ^ 


_  3(cQ,  )  _  3(eQ1 )  .  3(cQ1  ) 

'kz  '  S'CTT''  ~kz  3Q3'“  '  1  "”kz  35jp 


3(cQ1  ) 


3 (cQ.  ) 
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where 


,5  3(pJ)  s  3(CQl)  ^  3§><  v  ,  ,  ,  *  c 

5J  “  33] —  9k  “35] —  c<31  -3CJ*  J  ’  7* 2* 3> u’5- 

Let  denote  the  element  in  row  i  and  column  j  of  the  Jacobian 

and  let  (Kj_),  denote  element  i  of  the  split  flux  vector  .  Equation 
(3-3)  can  be  rewritten  as 


,  3(K  ) 

(Kl)ij  "  (Kl}i  3^3“  +  xk^“3;j- — )  * 


i  -  1.2, 3,4,5, 
J  -  1,2, 3, 4, 5, 
l  -  1,4,5. 


Explicitly  writing  all  the  elements  of  i?1  yields 
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(^1 ) 42 
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(Ki>54  "  -y-  tT  (j  kz  *  $7  Sk^ 

( K  ^ )  55  *  0 . 


Writing  out  all  25  elements  of  Rjj( 
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The  elements  of  R^  are 
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When  forming  the  Jacobian  of  or  Ka  in  Equation  (3-1),  the  deriv¬ 
atives  of  only  those  components  with  a  non-zero  coefficient  are  used. 
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